diff --git a/ShoulderCase/@Glenoid/calcDensity.m b/ShoulderCase/@Glenoid/calcDensity.m
old mode 100644
new mode 100755
index 8487779..3b6e81f
--- a/ShoulderCase/@Glenoid/calcDensity.m
+++ b/ShoulderCase/@Glenoid/calcDensity.m
@@ -1,903 +1,910 @@
 function outputArg = calcDensity(obj)
 %% Calculate the density of the glenoid bone in 6 volumes of interest
 
 % This function was developped by Paul Cimadomo during a semester  project (Fall 2018-19)
 % and implemented by Alexandre Terrier (2019-02-03)
 % sCase --> SCase
 % obj.cylinder removed
 % obj.volume removed
+% try-catch added for dicomreadVolume
 
 
-tic;
-disp('calcDensity');
+%tic;
+%disp('calcDensity');
+
 %Define a cylinder of 40mm height enclosing the glenoid
 
 SCase = obj.SCase;
 glenSurf = obj.surface.points;
 
 % Get local coordinate system from parent scapula object
 origin = SCase.shoulder.scapula.coordSys.origin;
 xAxis = SCase.shoulder.scapula.coordSys.PA;
 yAxis = SCase.shoulder.scapula.coordSys.IS;
 zAxis = SCase.shoulder.scapula.coordSys.ML;
 
 % Glenoid center is the closet point between the mean of the extracted
 % surface and all the points of the extracted surface.
 % This might be improved by find the "real" center, instead of
 % the closest one from the surface nodes.
 glenCenter = mean(glenSurf);
 vect = zeros(size(glenSurf));
 dist = zeros(size(glenSurf,1),1);
 for i=1:size(glenSurf,1)
     vect(i,:) = glenSurf(i,:) - glenCenter;
     dist(i) = norm(vect(i,:));
 end
 glenCenterNode = find(dist == min(dist));
 glenCenter = glenSurf(glenCenterNode,:);
 
 % Rotation matrix to align to scapula coordinate system
 R = [xAxis' yAxis', zAxis'];
 
 % Glenoid surface in scapular coordinate system
 glenSurf = (SCase.shoulder.scapula.glenoid.surface.points - origin) * R ;
 
 %Calculate the distance between the center of the scapular
 %coordinate system and each point of the glenoid surface
 
 glenoidCenterLocal = [obj.centerLocal.x,obj.centerLocal.y,obj.centerLocal.z];
 glenSurf = glenSurf - glenoidCenterLocal;
 
 %Create a unit vector in the direction of the cylinder axis
 unit_vector = zAxis;
 
 %Find the spinoglenoid notch located on the cylinder axis aligned with
 %its medial aspect of the cylinder with the
 spinoglenoid_notch = glenCenter + unit_vector * -glenoidCenterLocal(3) ;
 
 norm_vector=[];
 
 for n = 1:length(glenSurf);
     norm_vector = [norm_vector ; norm(glenSurf(n,1:2))];
 end
 
 %Define the radius of the cylinder as the maximum distance
 %between the scapular center and the glenoid points
 
 cylinder_radius = max(norm_vector) ;
 
 % obj.cylinder = [cylinder_radius, obj.center, obj.centerLine] ;
 sphereRadius = obj.radius ;
 sphereCenterGlen = obj.center + obj.centerLine ;
 %% Including cylinder in the images
 
 %Extract the images as a 4D matrix where spatial is the
 %positions of the pixels in the CT coordinate system
 dataCTPath = SCase.dataCTPath;
 filename = [dataCTPath '/dicom'];
-[V,spatial,dim] = dicomreadVolume(filename);
+try
+  [V,spatial,dim] = dicomreadVolume(filename);
+catch
+  outputArg = 1;
+  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
 list = dir(filename) ;
 size_list = size(list) ;
 dicominfopath = [filename '/' list(size_list(1)).name] ;
 dicom_information = dicominfo(dicominfopath);
 
 %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) ;
 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 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,:);
 x_max = origin_image(1)+PixelSpacings(1,1)*size_V(1);
 y_max = origin_image(2)+PixelSpacings(1,2)*size_V(2);
 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([origin_image(1) , x_max],[1 size_V(1)],1);
 a_i = coefficients_i(1);
 b_i = coefficients_i(2);
 
 coefficients_j = polyfit([origin_image(2) , y_max],[1 size_V(2)],1);
 a_j = coefficients_j(1);
 b_j = coefficients_j(2);
 
 coefficients_k = polyfit([origin_image(3) , z_max],[1 size_V(3)],1);
 a_k = coefficients_k(1);
 b_k = coefficients_k(2);
 
 %Pixel size : 0.31 mm x 0.31-0.51 mm x 0.51 mm
 dr = 0.31; %radius increment of the cylinder
 da = 0.25/cylinder_radius; %angle increment of the cylinder
 dz = 0.31; %height increment of the cylinder
 
 %Create two arrays : one empty (only 0) to fill with rescaling slope
 %and one with rescaling intercept to rescale value of the voxel
 %stored in matlab to its true HU value.
 V_cyl = zeros(size_V) ;
 V_rescale = zeros(size_V) ;
 
 h = 40; %height of the cylinder in [mm]
 
 for z = -5:dz:h ;     %Make a h + 5 mm cylinder to be able to
    for r = 0:dr:cylinder_radius ;
         for a = 0:da:2*pi ;
 
             x = r * cos(a) ;
             y = r * sin(a) ;
             z = z ;
 
             %Rotate and translate the cylinder to its position
             %in the CTs
             cylinder_coord_scap = R*[x;y;z] + spinoglenoid_notch';
 
             %Perform linear transformation CT (x,y,z) to image
             %(i,j,k) coordinate system
             i = a_i*cylinder_coord_scap(1) + b_i ;
             j = a_j*cylinder_coord_scap(2) + b_j ;
             k = a_k*cylinder_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
             V_cyl(j,i,k) = Rescale_slope ;
             V_rescale (j,i,k) = Rescale_intercept ;
             else
                 continue
             end
         end
     end
 end
 
 %Create a plate to close one side of the cylinder
 V_plate = zeros(size_V) ;
 for z = -5:0.1:-4.9;
    for y = -cylinder_radius:0.1:cylinder_radius ;
        for x = -cylinder_radius:0.1:cylinder_radius ;
 
             cylinder_coord_scap = R*[x;y;z] + spinoglenoid_notch';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*cylinder_coord_scap(1) + b_i ;
             j = a_j*cylinder_coord_scap(2) + b_j ;
             k = a_k*cylinder_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array of the plate with a value of
             %1000 to avoid erasing of the plate by thresholding
             V_plate(j,i,k) = 1000 ;
             else
                 continue
             end
        end
    end
 end
 
 V_plate_mask = ones(size_V) ;
 for z = -5:0.1:0;
    for y = -cylinder_radius:0.1:cylinder_radius ;
        for x = -cylinder_radius:0.1:cylinder_radius ;
 
             cylinder_coord_scap = R*[x;y;z] + spinoglenoid_notch';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*cylinder_coord_scap(1) + b_i ;
             j = a_j*cylinder_coord_scap(2) + b_j ;
             k = a_k*cylinder_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array of the plate with a value of
             %1000 to avoid erasing of the plate by thresholding
             V_plate_mask(j,i,k) = 0 ;
             else
                 continue
             end
 
        end
    end
 end
 
 %Create a sphere to exclude humeral head
 V_sphere = ones(size_V) ;
 for r = 0:dr:sphereRadius ;
    for a = 0:da/2:2*pi ;
        for phi = 0:da:pi;
 
             x = r * sin(phi) * cos(a) ;
             y = r * sin(phi) * sin(a) ;
             z = r * cos(phi);
 
             sphere_coord_scap = [x;y;z] + sphereCenterGlen';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*sphere_coord_scap(1) + b_i ;
             j = a_j*sphere_coord_scap(2) + b_j ;
             k = a_k*sphere_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
                 V_sphere(j,i,k) = 0;
             else
                 continue
             end
        end
    end
 end
 
 %Create masks of the 5 spherical shells
 V_sphere_SC = zeros(size_V) ;
 for r = 0:dr:sphereRadius + 3 ;
    for a = 0:da/2:2*pi ;
        for phi = 0:da:pi;
 
             x = r * sin(phi) * cos(a) ;
             y = r * sin(phi) * sin(a) ;
             z = r * cos(phi);
 
             sphere_coord_scap = [x;y;z] + sphereCenterGlen';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*sphere_coord_scap(1) + b_i ;
             j = a_j*sphere_coord_scap(2) + b_j ;
             k = a_k*sphere_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
             V_sphere_SC(j,i,k) = 1;
             else
                 continue
             end
        end
    end
 end
 
 V_sphere_ST = zeros(size_V) ;
 for r = 0:dr:sphereRadius + 6 ;
    for a = 0:da/2:2*pi ;
        for phi = 0:da:pi;
 
             x = r * sin(phi) * cos(a) ;
             y = r * sin(phi) * sin(a) ;
             z = r * cos(phi);
 
             sphere_coord_scap = [x;y;z] + sphereCenterGlen';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*sphere_coord_scap(1) + b_i ;
             j = a_j*sphere_coord_scap(2) + b_j ;
             k = a_k*sphere_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
             V_sphere_ST(j,i,k) = 1;
             else
                 continue
             end
        end
    end
 end
 
 V_sphere_T1 = zeros(size_V) ;
 for r = 0:dr:sphereRadius + 9;
    for a = 0:da/2:2*pi ;
        for phi = 0:da:pi;
 
             x = r * sin(phi) * cos(a) ;
             y = r * sin(phi) * sin(a) ;
             z = r * cos(phi);
 
             sphere_coord_scap = [x;y;z] + sphereCenterGlen';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*sphere_coord_scap(1) + b_i ;
             j = a_j*sphere_coord_scap(2) + b_j ;
             k = a_k*sphere_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
             V_sphere_T1(j,i,k) = 1;
             else
                 continue
             end
        end
    end
 end
 
 V_sphere_T2 = zeros(size_V) ;
 for r = 0:dr:sphereRadius + 12 ;
    for a = 0:da/2:2*pi ;
        for phi = 0:da:pi;
 
             x = r * sin(phi) * cos(a) ;
             y = r * sin(phi) * sin(a) ;
             z = r * cos(phi);
 
             sphere_coord_scap = [x;y;z] + sphereCenterGlen';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*sphere_coord_scap(1) + b_i ;
             j = a_j*sphere_coord_scap(2) + b_j ;
             k = a_k*sphere_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
             V_sphere_T2(j,i,k) = 1;
             else
                 continue
             end
        end
    end
 end
 
 V_sphere_T3 = zeros(size_V) ;
 for r = 0:dr:sphereRadius + 15;
    for a = 0:da/2:2*pi ;
        for phi = 0:da/3:pi;
 
             x = r * sin(phi) * cos(a) ;
             y = r * sin(phi) * sin(a) ;
             z = r * cos(phi);
 
             sphere_coord_scap = [x;y;z] + sphereCenterGlen';
 
             %Perform linear transformation CT to image
             %coordinate system
             i = a_i*sphere_coord_scap(1) + b_i ;
             j = a_j*sphere_coord_scap(2) + b_j ;
             k = a_k*sphere_coord_scap(3) + b_k ;
 
             %Round up to get integers
             i = round(i) ;
             j = round(j) ;
             k = round(k) ;
 
             %Ensure that the image coordinates don't exceed the
             %size of the original images (constrain it at the
             %bottom by avoiding negative cases and at the top
             %with the maximum number of pixels in one
             %direction)
             if (i > 0) && (j > 0) && (k > 0) && (i <= size_V(1)) && (j<=size_V(2)) && (k<=size_V(3))
             %Fill the empty array and mask with the position of
             %the cylinder
             V_sphere_T3(j,i,k) = 1;
             else
                 continue
             end
        end
    end
 end
 
 %Create an array filled with images + cylinder for
 %visualisation in plotSCase. 1000 multiplicator is chosen so
 %isosurface rendering can be used easily.
 %Get normal visualisation by setting z-axis to 1.5 units/voxel
 %(a_i/a_k = a_j/a_k) in volumeViewer
 
 V_visual = V+1000*V_cyl;  %can be used in volumeViewer for rendering
 
 %Create an array filled with the images contained within the
 %cylinder. Rescale it to the true HU values
 V_paint = V.*V_cyl + V_rescale;
 
 %Set to zero the domain covered by the sphere (so the humerus)
 V_paint = V_paint.*V_sphere;
 
 %Close the cylinder with the plate
 V_paint = V_paint + V_plate;
 
 
 %             %% Image processing
 %             %Suppose we take care of only one frame (frame 197 lets say).
 %
 %            I = V_paint(:,:,197) ;
 %            I_plate_mask = V_plate_mask(:,:,197) ;
 %
 %
 %            %Apply a gaussian filter to smooth the edges. Sigma indicates
 %            %the Gaussian smoothing kernels.
 %            sigma = 0.3;
 %            I_blurred = imgaussfilt(I,sigma) ;
 %
 %            %Set a threshold of 300 HU (take only bone)
 %            threshold = 300;
 %            I_threshold = I_blurred > threshold;
 %
 %            %Fill the holes in the image with 8-connected pixel for filling
 %            I_filled = imfill(I_threshold,8,'holes');
 %
 %            %Clean the remaining parts that are outside the glenoid part we consider
 %            I_cleaned = bwareaopen(I_filled,50,4) ;
 %
 %            %Create the image where the glenoid bone is filled and substract
 %            %from it the binarized plate to not take it into account in the
 %            %calculation of the mean
 %            I_plate_bin = imbinarize(V_plate(:,:,197));
 %            I_mask = imbinarize(I_cleaned - I_plate_bin) ;
 %
 %            %Create a structural element having a disk shape with a radius
 %            %of the wanted erosion. Then convert it to image coordinate system
 %            %by dividing by the spacing in mm between two pixels
 %            erosion = 3 ;%[mm]
 %            erosion_image_coord = floor(erosion/PixelSpacings(1,1)) ;
 %            SE = strel('disk', erosion_image_coord, 8) ; %
 %            I_eroded = imerode(I_mask, SE) ;
 %
 %            %Separate the image between cortical bone and trabecular
 %            I_mask_fin = bwareaopen(I_mask.*I_plate_mask,50,4);
 %            I_eroded_fin = I_eroded.*I_plate_mask ;
 %
 %            I_cortical = (I_mask_fin-I_eroded_fin) ;
 %
 %            %Get SC, ST, T1, T2 and T3: start with SC sphere and then
 %            %substract in the right order the different VOIs
 %
 %            I_SC = imfill(V_sphere_SC(:,:,197),8,'holes') ;
 %            I_ST = imfill(V_sphere_ST(:,:,197),8,'holes') - imfill(V_sphere_SC(:,:,197),8,'holes');
 %            I_T1 = imfill(V_sphere_T1(:,:,197),8,'holes') - imfill(V_sphere_ST(:,:,197),8,'holes') ;
 %            I_T2 = imfill(V_sphere_T2(:,:,197),8,'holes') - imfill(V_sphere_T1(:,:,197),8,'holes') ;
 %            I_T3 = imfill(V_sphere_T3(:,:,197),8,'holes') - imfill(V_sphere_T2(:,:,197),8,'holes') ;
 %
 %            I_SC_fin = imbinarize(I_SC.*I_mask_fin) ;
 %            I_ST_fin = imbinarize(I_ST.*I_eroded_fin) ;
 %            I_T1_fin = imbinarize(I_T1.*I_eroded_fin) ;
 %            I_T2_fin = imbinarize(I_T2.*I_eroded_fin) ;
 %            I_T3_fin = imbinarize(I_T3.*I_eroded_fin) ;
 %            I_CO_fin = imbinarize(I_cortical - I_SC_fin) ;
 %
 % %            %Flatten the image to ease the calculation of the mean
 %            size_I = size(I) ;
 %            size_mask = size(I_mask) ;
 % %
 %            I_mask_flat = reshape(I_mask_fin,size_mask(1)*size_mask(2),1) ;
 %            I_flat = reshape(I,size_I(1)*size_I(2),1) ;
 %
 %            mean_vector = [];
 %            count_HU = 0;
 %
 %            %For value in the flat mask = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_mask_flat(i) == 1
 %                   mean_vector = [mean_vector ; I_flat(i)] ;
 %                   count_HU = count_HU + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean = sum(mean_vector)/count_HU ;
 %            %% Measure mean in the CO
 %
 %            I_CO_flat = reshape(I_CO_fin ,size_mask(1)*size_mask(2),1) ;
 %
 %            mean_vector_CO = [];
 %            count_HU_CO = 0;
 %
 %            %For value in the flat SC = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_CO_flat(i) == 1
 %                   mean_vector_CO = [mean_vector_CO ; I_flat(i)] ;
 %                   count_HU_CO = count_HU_CO + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean_CO = sum(mean_vector_CO)/count_HU_CO ;
 %            %% Measure mean in the SC
 %
 %            I_SC_flat = reshape(I_SC_fin ,size_mask(1)*size_mask(2),1) ;
 %
 %            mean_vector_SC = [];
 %            count_HU_SC = 0;
 %
 %            %For value in the flat SC = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_SC_flat(i) == 1
 %                   mean_vector_SC = [mean_vector_SC ; I_flat(i)] ;
 %                   count_HU_SC = count_HU_SC + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean_SC = sum(mean_vector_SC)/count_HU_SC ;
 %
 %            %% Measure mean in the ST
 %
 %            I_ST_flat = reshape(I_ST_fin ,size_mask(1)*size_mask(2),1) ;
 %
 %            mean_vector_ST = [];
 %            count_HU_ST = 0;
 %
 %            %For value in the flat ST = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_ST_flat(i) == 1
 %                   mean_vector_ST = [mean_vector_ST ; I_flat(i)] ;
 %                   count_HU_ST = count_HU_ST + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean_ST = sum(mean_vector_ST)/count_HU_ST ;
 %
 %            %% Measure mean in the T1
 %            I_T1_flat = reshape(I_T1_fin ,size_mask(1)*size_mask(2),1) ;
 %
 %            mean_vector_T1 = [];
 %            count_HU_T1 = 0;
 %
 %            %For value in the flat T1 = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_T1_flat(i) == 1
 %                   mean_vector_T1 = [mean_vector_T1 ; I_flat(i)] ;
 %                   count_HU_T1 = count_HU_T1 + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean_T1 = sum(mean_vector_T1)/count_HU_T1 ;
 %
 %            %% Measure mean in the T2
 %            I_T2_flat = reshape(I_T2_fin ,size_mask(1)*size_mask(2),1) ;
 %
 %            mean_vector_T2 = [];
 %            count_HU_T2 = 0;
 %
 %            %For value in the flat T2 = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_T2_flat(i) == 1
 %                   mean_vector_T2 = [mean_vector_T2 ; I_flat(i)] ;
 %                   count_HU_T2 = count_HU_T2 + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean_T2 = sum(mean_vector_T2)/count_HU_T2 ;
 %
 %            %% Measure mean in the T3
 %            I_T3_flat = reshape(I_T3_fin ,size_mask(1)*size_mask(2),1) ;
 %
 %            mean_vector_T3 = [];
 %            count_HU_T3 = 0;
 %
 %            %For value in the flat T3 = 1 (where there is glenoid), set
 %            %the mean_vector to the value of the corresponding pixel value
 %            %of the glenoid image
 %            for i = 1:1:size_mask(1)*size_mask(2) ;
 %
 %                if I_T3_flat(i) == 1
 %                   mean_vector_T3 = [mean_vector_T3 ; I_flat(i)] ;
 %                   count_HU_T3 = count_HU_T3 + 1 ;
 %                else
 %                    continue
 %                end
 %            end
 %
 %            image_mean_T3 = sum(mean_vector_T3)/count_HU_T3 ;
 
 %% Volume processing
 %Initialize the mean vectors to calculate the mean in each part
 %of interest
 
 mean_vector_CO = [] ;
 mean_vector_SC = [];
 mean_vector_ST = [];
 mean_vector_T1 = [];
 mean_vector_T2 = [];
 mean_vector_T3 = [];
 
 % For loop across all the images to calculate glenoid bone
 % quality in each VOI of each image
 for n = 1:1:size_V(3) ;
 
 %Focus on one image and corresponding mask
 I = V_paint(:,:,n) ;
 I_plate_mask = V_plate_mask(:,:,n) ;
 
 %Apply a gaussian filter to smooth the edges. Sigma indicates
 %the Gaussian smoothing kernels.
 sigma = 0.3;
 I_blurred = imgaussfilt(I,sigma) ;
 
 %Segment bone with a threshold of 300 HU
 threshold = 300;
 I_threshold = I_blurred > threshold;
 
 %Fill the holes in the image with 8-connected pixel for filling
 I_filled = imfill(I_threshold,8,'holes');
 
 %Clean the remaining parts that are outside the glenoid part we consider
 I_cleaned = bwareaopen(I_filled,50,4) ;
 
 %Create the image where the glenoid bone is filled and substract
 %from it the binarized plate to not take it into account in the
 %calculation of the mean
 I_plate_bin = imbinarize(V_plate(:,:,n));
 I_mask = imbinarize(I_cleaned - I_plate_bin) ;
 
 %Create a structural element having a disk shape with a radius
 %of the wanted erosion. Then convert it to image coordinate system
 %by dividing by the spacing in mm between two pixels
 erosion = 3 ;%[mm]
 erosion_image_coord = floor(erosion/PixelSpacings(1,1)) ;
 SE = strel('disk', erosion_image_coord, 8) ; %
 I_eroded = imerode(I_mask, SE) ;
 
 %Separate the image between cortical bone and trabecular
 I_mask_fin = bwareaopen(I_mask.*I_plate_mask,50,4);
 I_eroded_fin = I_eroded.*I_plate_mask ;
 
 I_cortical = (I_mask_fin-I_eroded_fin) ;
 
 %Get SC, ST, T1, T2 and T3 by substracting each VOI by the
 %previous one
 
 I_SC = imfill(V_sphere_SC(:,:,n),8,'holes') ;
 I_ST = imfill(V_sphere_ST(:,:,n),8,'holes') - imfill(V_sphere_SC(:,:,n),8,'holes');
 I_T1 = imfill(V_sphere_T1(:,:,n),8,'holes') - imfill(V_sphere_ST(:,:,n),8,'holes') ;
 I_T2 = imfill(V_sphere_T2(:,:,n),8,'holes') - imfill(V_sphere_T1(:,:,n),8,'holes') ;
 I_T3 = imfill(V_sphere_T3(:,:,n),8,'holes') - imfill(V_sphere_T2(:,:,n),8,'holes') ;
 
 %Binarize the image to only have 0 and 1 in the final VOI masks
 I_SC_fin = imbinarize(I_SC.*I_mask_fin) ;
 I_ST_fin = imbinarize(I_ST.*I_eroded_fin) ;
 I_T1_fin = imbinarize(I_T1.*I_eroded_fin) ;
 I_T2_fin = imbinarize(I_T2.*I_eroded_fin) ;
 I_T3_fin = imbinarize(I_T3.*I_eroded_fin) ;
 I_CO_fin = imbinarize(I_cortical - I_SC_fin) ;
 
 %Flatten the image to ease the calculation of the mean
 size_I = size(I) ;
 size_mask = size(I_mask) ;
 
 I_flat = reshape(I,size_I(1)*size_I(2),1) ;
 
 
 %% Measure mean in the CO
 
 I_CO_flat = reshape(I_CO_fin ,size_mask(1)*size_mask(2),1) ;
 
 
 %For value in the flat CO = 1 (where there is glenoid), set
 %the mean_vector to the value of the corresponding pixel value
 %of the glenoid image
 for i = 1:1:size_mask(1)*size_mask(2) ;
 
    if I_CO_flat(i) == 1
       mean_vector_CO = [mean_vector_CO ; I_flat(i)] ;
    else
        continue
    end
 end
 
 %% Measure mean in the SC
 
 I_SC_flat = reshape(I_SC_fin ,size_mask(1)*size_mask(2),1) ;
 
 
 %For value in the flat SC = 1 (where there is glenoid), set
 %the mean_vector to the value of the corresponding pixel value
 %of the glenoid image
 for i = 1:1:size_mask(1)*size_mask(2) ;
 
    if I_SC_flat(i) == 1
       mean_vector_SC = [mean_vector_SC ; I_flat(i)] ;
    else
        continue
    end
 end
 
 
 %% Measure mean in the ST
 
 I_ST_flat = reshape(I_ST_fin ,size_mask(1)*size_mask(2),1) ;
 
 %For value in the flat ST = 1 (where there is glenoid), set
 %the mean_vector to the value of the corresponding pixel value
 %of the glenoid image
 for i = 1:1:size_mask(1)*size_mask(2) ;
 
    if I_ST_flat(i) == 1
       mean_vector_ST = [mean_vector_ST ; I_flat(i)] ;
    else
        continue
    end
 end
 
 %% Measure mean in the T1
 I_T1_flat = reshape(I_T1_fin ,size_mask(1)*size_mask(2),1) ;
 
 %For value in the flat T1 = 1 (where there is glenoid), set
 %the mean_vector to the value of the corresponding pixel value
 %of the glenoid image
 for i = 1:1:size_mask(1)*size_mask(2) ;
 
    if I_T1_flat(i) == 1
       mean_vector_T1 = [mean_vector_T1 ; I_flat(i)] ;
    else
        continue
    end
 end
 
 %% Measure mean in the T2
 I_T2_flat = reshape(I_T2_fin ,size_mask(1)*size_mask(2),1) ;
 
 %For value in the flat T2 = 1 (where there is glenoid), set
 %the mean_vector to the value of the corresponding pixel value
 %of the glenoid image
 for i = 1:1:size_mask(1)*size_mask(2) ;
 
    if I_T2_flat(i) == 1
       mean_vector_T2 = [mean_vector_T2 ; I_flat(i)] ;
    else
        continue
    end
 end
 
 %% Measure mean in the T3
 I_T3_flat = reshape(I_T3_fin ,size_mask(1)*size_mask(2),1) ;
 
 %For value in the flat T3 = 1 (where there is glenoid), set
 %the mean_vector to the value of the corresponding pixel value
 %of the glenoid image
 for i = 1:1:size_mask(1)*size_mask(2) ;
 
    if I_T3_flat(i) == 1
       mean_vector_T3 = [mean_vector_T3 ; I_flat(i)] ;
    else
        continue
    end
 end
 
 end
 
-mean_CO = mean(mean_vector_CO)
-mean_SC = mean(mean_vector_SC)
-mean_ST = mean(mean_vector_ST)
-mean_T1 = mean(mean_vector_T1)
-mean_T2 = mean(mean_vector_T2)
-mean_T3 = mean(mean_vector_T3)
+mean_CO = mean(mean_vector_CO);
+mean_SC = mean(mean_vector_SC);
+mean_ST = mean(mean_vector_ST);
+mean_T1 = mean(mean_vector_T1);
+mean_T2 = mean(mean_vector_T2);
+mean_T3 = mean(mean_vector_T3);
 
 %obj.volume = V_visual;
 
 obj.density = [mean_CO, mean_SC, mean_ST, mean_T1, mean_T2, mean_T3];
-toc
+%toc
 
 outputArg = 1;
 end
diff --git a/measureSCase.m b/measureSCase.m
index 14bc9b0..6cf04df 100755
--- a/measureSCase.m
+++ b/measureSCase.m
@@ -1,450 +1,452 @@
 function [status,message] = measureSCase(varargin)
 % MEASURESCASE Update the entire SCase DB with anatimical measurements.
 
 % The function can be run from (lbovenus) server with the following shell command
 % cd /home/shoulder/methods/matlab/database; matlab -nodisplay -nosplash -r "measureSCase;quit"
 
 % Progress can be checked through log file with following shell command
 % cd /home/shoulder/methods/matlab/database;tail -f log/measureSCase.log
 
 % It take ~30 min for 700 cases when run from server.
 % After run from server, the output file
 % /shoulder/data/Excel/xlsFromMatlab/anatomy.csv
 % must be open from Excel and saved
 % as xls at the same location.
 
 % The script measureSCase will replace 3 previous functions
 % 1) rebuildDatabaseScapulaMeasurements.m --> list of all SCases --> execute scapula_measure.m
 % 2) scapula_measure.m --> execute scapula_calculation and save results in Excel and SQL
 % 3) scapula_calculation.m --> perform anatomical analysis --> results varaiable
 
 % The current version produces correctly the entire csv when run
 % without arguments only (for the entire database)
 
 % Input
 %    Either empty, 'N', 'P', or a SCase ('P315')
 
 % Output
 %    File /shoulder/data/Excel/xlsFromMatlab/anatomy.csv
 %    File /shoulder/data/Excel/xlsFromMatlab/anatomy.xls (only windows)
 %    File /home/shoulder/data/matlab/SCaseDB.mat
 %    File {SCaseId}.mat in each {SCaseId}/{CT}/matlab directory
 %    logs in /log/measureSCase.log
 
 % Example: measureSCase
 
 % Author: Alexandre Terrier, EPFL-LBO
 % Creation date: 2018-07-01
 % Revision date: 2018-12-31
 
 % TO DO:
 %  Read first last anatomy.csv and/or SCaseDB.mat
 %  Fill with patient and clinical data, from excel or MySQL
 %  use argument to redo or update
 %  Save csv within the main loop (not sure it's a good idea)
 %  Add more message in log
 %  Add more test to assess validity of data
 %  Update MySQL
 %  Rename all objects with a capital
 %  Make it faster by not loading scapula surface if already in SCaseDB.mat
 %  Add argument 'load' to force reload files even if they are alerady in DB
 %  Add argument 'measure' to force re-measuring event if measurement is already in DB
 %  Add "scapula_addmeasure.m" --> should be coded in ShoulderCase class
 
 
 tic; % Start stopwatch timer
 
 %% Open log file
 logFileID = openLogFile('measureSCase.log');
 
 %% Set the data directory from the configuration file config.txt
 dataDir = openConfigFile('config.txt', logFileID);
 
 
 
 % configFile = 'config.txt';
 % if ~exist(configFile, 'file')
 %     error([configFile, ' required']);
 % end
 % fileID = fopen(configFile,'r');
 % dataDir = fscanf(fileID, '%s'); % Read config file without spaces
 % fclose(fileID);
 % k = strfind(dataDir, 'dataDir='); % Find definition of dataDir
 % k = k + 8; % Remove 'dataDir='
 % dataDir = dataDir(k:end); % Assume that config file ends with dataDir content
 % % Check if dataDir exists
 % if ~exist(dataDir, 'dir')
 %     fprintf(logFileID, ['Data directory not found, check ' configFile]);
 %     error(['Data directory not found, check ' configFile]);
 % end
 
 
 %% Location of the XLS ShoulderDatabase
 xlsDir    = [dataDir '/Excel/xlsFromMatlab'];
 matlabDir = [dataDir '/matlab'];
 
 %% Get list of all SCase
 % Get list of SCases from varargin
 fprintf(logFileID, '\n\nList of SCase');
 if nargin == 0 % findScapulaLandmarksSCase called without argument
     SCaseList = listSCase;
 elseif nargin == 1
     inputArg = varargin{1};
     SCaseList = listSCase(inputArg);
 else
     error('inputArg can be empty, N, P, or a SCaseID')
 end
 
 %% Load current xls database (for clinical data)
 fprintf(logFileID, '\nLoad xls database');
 addpath('XLS_MySQL_DB');
 filename = [dataDir '/Excel/ShoulderDataBase.xlsx'];
 excelRaw = rawFromExcel(filename); % cell array
 
 excelSCaseID = excelRaw(2:end, 1);
 excelDiagnosis = excelRaw(2:end, 4);
 excelPatientHeight = excelRaw(2:end, 23);
 excelGlenoidWalch = excelRaw(2:end, 55);
 
 % Lines below adapted from XLS_MySQL_DB/MainExcel2XLS
 % Excel.diagnosisList = diagnosisListFromExcel(Excel.Raw);
 % Excel.treatmentList = treatmentListFromExcel(Excel.Raw);
 % ExcelData.Patient = patientFromExcel(excelRaw); % Structure with patient data
 % Excel.shoulder = shoulderFromExcel(Excel.patient, Excel.Raw);
 % [Excel.SCase, Excel.diagnosis, Excel.treatment, Excel.outcome, Excel.study] = SCaseFromExcel(...
 %     Excel.shoulder, Excel.patient, Excel.diagnosisList, Excel.treatmentList, Excel.Raw);
 % fprintf(logFileID, ': OK');
 
 %% Add path to ShoulderCase class
 addpath('ShoulderCase');
 
 %% Instance of a ShoulderCase object
 if (exist('SCase','var') > 0)
     clear SCase; % Check for delete method
 end
 SCase = ShoulderCase; % Instanciate a ShoulderCase object
 SCase.dataPath = dataDir; % Set dataDir for SCase
 
 %% Start loop over SCases in SCaseList
 nSCase = length(SCaseList); % Number of SCases
 for iSCaseID = 1:nSCase
     SCaseID = SCaseList(iSCaseID).id; %
     SCase(iSCaseID).id = SCaseID;
     SCase(iSCaseID).dataPath = dataDir; % Set dataDir for SCase
     percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f');
     fprintf(logFileID, ['\n\nSCaseID: ' SCaseID ' (' percentProgress '%%)']);
 
     % There are 3 parts within this SCase loop:
     % 1) Load & analyses manual data from amira
     % 2) Load & analyses auto   data from matlab (Statistical Shape Model)
     % 3) Load clinical data from Excel database, and set them to SCase
     % 4) Save SCase in file SCaseID/matlab/SCase.mat
 
     % 1) Load & analyses manual data from amira
     fprintf(logFileID, '\n  Segmentation manual ');
     output = SCase(iSCaseID).path; % Set data path of this SCase
     if output % Continue if amira dir exists in SCase
         fprintf(logFileID, '\n    Load scapula surface and landmarks');
         output = SCase(iSCaseID).shoulder.scapula.load;
         if output
             fprintf(logFileID, ': OK');
             fprintf(logFileID, '\n    Set coord. syst.');
             output = SCase(iSCaseID).shoulder.scapula.coordSysSet;
             if output
                 fprintf(logFileID, ': OK');
                 fprintf(logFileID, '\n    Load glenoid surface');
                 output = SCase(iSCaseID).shoulder.scapula.glenoid.load;
                 if output
                     fprintf(logFileID, ': OK');
                     fprintf(logFileID, '\n    Measure glenoid anatomy');
                     output = SCase(iSCaseID).shoulder.scapula.glenoid.morphology;
                     if output
                         fprintf(logFileID, ': OK');
                         fprintf(logFileID, '\n    Measure glenoid density');
                         output = SCase(iSCaseID).shoulder.scapula.glenoid.calcDensity;
                         if output
                           fprintf(logFileID, ': OK');
+                        else
+                          fprintf(logFileID, ': Could no read dicom');
                         end
                         fprintf(logFileID, '\n    Load humerus data');
                         output = SCase(iSCaseID).shoulder.humerus.load;
                         if output
                             fprintf(logFileID, ': OK');
                             fprintf(logFileID, '\n    Measure humerus subluxation');
                             output = SCase(iSCaseID).shoulder.humerus.subluxation;
                             if output
                                 fprintf(logFileID, ': OK');
                                 fprintf(logFileID, '\n    Measure acromion anatomy');
                                 % Should be run after SCase.shoulder.humerus (for AI)
                                 output = SCase(iSCaseID).shoulder.scapula.acromion.morphology;
                                 if output
                                     fprintf(logFileID, ': OK');
                                 else
                                     fprintf(logFileID, ': Error');
                                 end % Acromion anatmy
                             else
                                 fprintf(logFileID, ': Error');
                             end % Subluxation humerus
                         else
                             fprintf(logFileID, ': Error');
                         end % Load humerus
                     else
                         fprintf(logFileID, ': Error');
                     end  % Anatomy glenoid
                 else
                     fprintf(logFileID, ': Error');
                 end  % Load glenoid
             else
                 fprintf(logFileID, '\n    No amira folder');
             end % Set coord. syst.
         end % Load scapula
     end % Amira path exist
 
     % 2) Load & analyses auto data from matlab (Statistical Shape Model)
     fprintf(logFileID, '\n  Segmentation auto ');
     % Load scapula surface and landmarks
     fprintf(logFileID, '\n    Load scapula surface and landmarks');
     output = SCase(iSCaseID).shoulder.scapulaAuto.loadAuto; % Load auto data from matlab directory
     if output
         fprintf(logFileID, ': OK');
         fprintf(logFileID, '\n    Set coord. syst.');
         output = SCase(iSCaseID).shoulder.scapulaAuto.coordSysSet;
         if output
             fprintf(logFileID, ': OK');
             fprintf(logFileID, '\n    Load glenoid surface');
             output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.loadAuto;
             if output
                 fprintf(logFileID, ': OK');
                 fprintf(logFileID, '\n    Measure glenoid anatomy');
                 output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.morphology;
                 if output
                     fprintf(logFileID, ': OK');
                 else
                     fprintf(logFileID, ': Error');
                 end
             else
                 fprintf(logFileID, ': Error');
             end
         else
             fprintf(logFileID, ': Error');
         end
     else
         fprintf(logFileID, ': Error');
     end
 
     % 3) Set clinical data from Excel database to SCase
     fprintf(logFileID, '\n  Set clinical data from Excel');
 
     % Get idx of excelRaw for SCaseID
     % Use cell2struct when dots are removed from excel headers
     idx = find(strcmp(excelRaw, SCaseID)); % Index of SCaseID in excelRaw
     SCase(iSCaseID).diagnosis = excelRaw{idx,4};
     SCase(iSCaseID).treatment = excelRaw{idx,9};
     SCase(iSCaseID).patient.gender = excelRaw{idx,19};
     SCase(iSCaseID).patient.age = excelRaw{idx,21};
     SCase(iSCaseID).patient.height = excelRaw{idx,23};
     SCase(iSCaseID).patient.weight = excelRaw{idx,24};
     SCase(iSCaseID).patient.BMI = excelRaw{idx,25};
     SCase(iSCaseID).shoulder.scapula.glenoid.walch = excelRaw{idx,55};
 
 
     % Patient data --> TO DO
     % find patient id (in SQL) corresponding to iSCaseID
 %
 %
 %     idxSCase = find(contains(ExcelData, SCaseID));
 %     SCase(iSCaseID).patient.id        = ExcelData(idxSCase).patient.patient_id; % SQL patient id
 %     SCase(iSCaseID).patient.idMed     = ExcelData(idxSCase).patient.IPP; % IPP to be replace by coded (anonymized) IPP from SecuTrial
 %     SCase(iSCaseID).patient.gender    = ExcelData(idxSCase).patient.gender;
 %     SCase(iSCaseID).patient.age       = []; % Age (years) at treatment, or preop CT (we might also add birthdate wi day et at 15 of the month)
 %     SCase(iSCaseID).patient.ethnicity = ExcelData(idxSCase).patient.IPP;
 %     SCase(iSCaseID).patient.weight    = ExcelData(idxSCase).patient.weight;
 %     SCase(iSCaseID).patient.height    = ExcelData(idxSCase).patient.height;
 %     SCase(iSCaseID).patient.BMI       = ExcelData(idxSCase).pateint.BMI;
 %     SCase(iSCaseID).patient.comment   = ExcelData(idxSCase).patient.comment;
 
     % 4) Save SCase in file SCaseID/matlab/SCase.mat
     fprintf(logFileID, '\n  Save SCase');
     output = SCase(iSCaseID).saveMatlab;
     if output
         fprintf(logFileID, ': OK');
     else
         fprintf(logFileID, ': Error');
     end
 end % End of loop on SCaseList
 
 %% Save the entire SCase array as a matlab file
 fprintf(logFileID, '\n\n  Save SCase database');
 filename = 'SCaseDB';
 filename = [matlabDir '/' filename '.mat'];
 try
     SCaseDB = SCase;
     save(filename, 'SCaseDB');
     fprintf(logFileID, ': OK');
 catch
     fprintf(logFileID, ': Error');
 end
 
 %% Write csv file
 % This might be a function (SCase.csvSave()). The input would be a filename and a structure
 % data
 
 fprintf(logFileID, '\n  Save csv database');
 
 % Replace header and data by structure. Currently not working
 %{
 txtFilename = [xlsDir, '/anatomy.txt']; % Name of the txt file
 
 DataStruc = struc(...
     'SCase_id',              SCase.id, ...
     'shoulder_side',         SCase.shoulder.side,...
     'glenoid_radius',        SCase.shoulder.scapula.glenoid.radius,...
     'glenoid_sphereRMSE',    SCase.shoulder.scapula.glenoid.sphereRMSE,...
     'glenoid_depth',         SCase.shoulder.scapula.glenoid.depth,...
     'glenoid_width',         SCase.shoulder.scapula.glenoid.width,...
     'glenoid_height',        SCase.shoulder.scapula.glenoid.height,...
     'glenoid_centerPA',      SCase.shoulder.scapula.glenoid.centerLocal(1),...
     'glenoid_centerIS',      SCase.shoulder.scapula.glenoid.centerLocal(2),...
     'glenoid_centerML',      SCase.shoulder.scapula.glenoid.centerLocal(3),...
     'glenoid_versionAmpl',   SCase.shoulder.scapula.glenoid.versionAmpl,...
     'glenoid_versionOrient', SCase.shoulder.scapula.glenoid.versionOrient,...
     'glenoid_version',       SCase.shoulder.scapula.glenoid.version,...
     'glenoid_inclination',   SCase.shoulder.scapula.glenoid.inclination,...
     'humerus_jointRadius',   SCase.shoulder.humerus.jointRadius,...
     'humerus_headRadius',    SCase.shoulder.humerus.radius,...
     'humerus_GHSAmpl',       SCase.shoulder.humerus.GHSAmpl,...
     'humerus_GHSOrient',     SCase.shoulder.humerus.GHSOrient,...
     'humerus_SHSAmpl',       SCase.shoulder.humerus.SHSAmpl,...
     'humerus_SHSOrient',     SCase.shoulder.humerus.SHSOrient,...
     'humerus_SHSAngle',      SCase.shoulder.humerus.SHSAngle,...
     'humerus_SHSPA',         SCase.shoulder.humerus.SHSPA,...
     'humerus_SHSIS',         SCase.shoulder.humerus.SHSIS,...
     'acromion_AI',           SCase.shoulder.scapula.acromion.AI,...
     'acromion_CSA',          SCase.shoulder.scapula.acromion.CSA,...
     'acromion_PS',           SCase.shoulder.scapula.acromion.PS...
     );
 DataTable = strct2table(Data);
 writetable(DataTable,filename);
 %}
 
 % Header of the csv file
 dataHeader = [...
     'SCase_id,' ...
     'shoulder_side,' ...
     'glenoid_radius,' ...
     'glenoid_sphereRMSE,' ...
     'glenoid_depth,' ...
     'glenoid_width,' ...
     'glenoid_height,' ...
     'glenoid_centerPA,' ...
     'glenoid_centerIS,' ...
     'glenoid_centerML,' ...
     'glenoid_versionAmpl,' ...
     'glenoid_versionOrient,' ...
     'glenoid_version,' ...
     'glenoid_inclination,' ...
     'humerus_jointRadius,' ...
     'humerus_headRadius,' ...
     'humerus_GHSAmpl,' ...
     'humerus_GHSOrient,' ...
     'humerus_SHSAmpl,' ...
     'humerus_SHSOrient,' ...
     'humerus_SHSAngle,' ...
     'humerus_SHSPA,' ...
     'humerus_SHSIS,' ...
     'acromion_AI,' ...
     'acromion_CSA,' ...
     'acromion_PSA,'...
     'acromion_AAA'...
     ];
 
 csvFilename = [xlsDir, '/anatomy.csv']; % Name of the csv file
 csvFileId = fopen(csvFilename, 'w');
 fprintf(csvFileId,dataHeader); % Write header
 for iSCaseID = 1:length(SCaseList) % Loop on SCaseList to write data in csv
     disp(SCase(iSCaseID).id); % For debug
     fprintf(csvFileId,['\n'...
         SCase(iSCaseID).id ','...                                               % SCase_id
         SCase(iSCaseID).shoulder.side ','...                                    % shoulder_side
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.radius) ','...         % glenoid_radius
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.sphereRMSE) ','...     % glenoid_sphereRMSE
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.depth) ','...          % glenoid_depth
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.width) ','...          % glenoid_width
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.height) ','...         % glenoid_height
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.centerLocal.x) ','... % glenoid_centerPA
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.centerLocal.y) ','... % glenoid_centerIS
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.centerLocal.z) ','... % glenoid_centerML
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.versionAmpl) ','...    % glenoid_versionAmpl
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.versionOrient) ','...  % glenoid_versionOrient
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.version) ','...        % glenoid_version
         num2str(SCase(iSCaseID).shoulder.scapula.glenoid.inclination) ','...    % glenoid_inclination
         num2str(SCase(iSCaseID).shoulder.humerus.jointRadius) ','...            % humerus_jointRadius
         num2str(SCase(iSCaseID).shoulder.humerus.radius) ','...                 % humerus_headRadius
         num2str(SCase(iSCaseID).shoulder.humerus.GHSAmpl) ','...                % humerus_GHSAmpl
         num2str(SCase(iSCaseID).shoulder.humerus.GHSOrient) ','...              % humerus_GHSOrient
         num2str(SCase(iSCaseID).shoulder.humerus.SHSAmpl) ','...                % humerus_SHSAmpl
         num2str(SCase(iSCaseID).shoulder.humerus.SHSOrient) ','...              % humerus_SHSOrient
         num2str(SCase(iSCaseID).shoulder.humerus.SHSAngle) ','...               % humerus_SHSAgle
         num2str(SCase(iSCaseID).shoulder.humerus.SHSPA) ','...                  % humerus_SHSPA
         num2str(SCase(iSCaseID).shoulder.humerus.SHSIS) ','...                  % humerus_SHSIS
         num2str(SCase(iSCaseID).shoulder.scapula.acromion.AI) ','...            % acromion_AI
         num2str(SCase(iSCaseID).shoulder.scapula.acromion.CSA) ','...           % acromion_CSA
         num2str(SCase(iSCaseID).shoulder.scapula.acromion.PSA) ','...            % acromion_PSA
         num2str(SCase(iSCaseID).shoulder.scapula.acromion.AAA)...                % acromion_AAA
         ]);
 end
 fclose(csvFileId);
 
 fprintf(logFileID, [': ' csvFilename]);
 
 %% Write xls file  (Windows only)
 [~,message] = csv2xlsSCase(csvFilename);
 fprintf(logFileID, message);
 
 % If run from non-windows system, only csv will be produced. The xls can be
 % produced by opening the csv from Excel and save as xls.
 % Could be a function with 1 input (cvsFilenames
 
 
 
 %{
 xlsFilename = strrep(csvFilename, '.csv', '.xls'); % Replace .csv by .xls
 if ispc % Check if run from windows!
     % Read cvs file as table and write the table as xls
     cvsTable = readtable(csvFilename); % Get table from csv (only way to read non-numeric data)
     cvsCell = table2cell(cvsTable); % Tranform table cell array
     dataHeader = cvsTable.Properties.VariableNames; % Get header
     cvsCell = [dataHeader; cvsCell]; % Add header
 
     sheet = 'anatomy'; % Sheet name of the xls file
 
     [status,message] = xlswrite(xlsFilename, cvsCell, sheet); % Save xls
     if status
         fprintf(logFileId, ['\nSaved ' xlsFilename]);
     else
         fprintf(logFileId, ['\nCould not save ' xlsFilename]);
         fprintf(logFileId, ['\n' message.identifier]);
         fprintf(logFileId, ['\n' message.message]);
     end
 else
     fprintf(logFileId, ['\n' xlsFilename ' not saved. Open and save as xls from Excel']);
 end
 %}
 
 
 %% Close log file
 fprintf(logFileID, '\n\nSCase measured: %s', num2str(length(SCaseList))); % Number of SCases measured
 fprintf(logFileID, '\nElapsed time (s): %s', num2str(toc)); % stopwatch timer
 fprintf(logFileID, '\n');
 fclose(logFileID); % Close log file
 
 %% Output of the SCase if required by argument
 % SCase.output % inputArg = abaqus/density/references/
 
 % update mySQL
 % SCase.sqlSave
 
 
 status  = 1;
 message = 'OK';
 
 end