diff --git a/ShoulderCase/@Glenoid/calcDensity.m b/ShoulderCase/@Glenoid/calcDensity.m
index 9a4bc62..ab9595e 100755
--- a/ShoulderCase/@Glenoid/calcDensity.m
+++ b/ShoulderCase/@Glenoid/calcDensity.m
@@ -1,910 +1,916 @@
 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');
 
 %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'];
 try
   [V,spatial,dim] = dicomreadVolume(filename);
 catch
   outputArg = false;
   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);
+dicominfopath = [filename '/' list(round(size_list(1)/2)).name] ;
+
+try                                                                             %
+  dicom_information = dicominfo(dicominfopath);                                 %   A try/catch for dicominfo()
+catch                                                                           %   has been added in case
+  outputArg = false;                                                            %   dicominfopath is not valid
+  return;                                                                       %
+end                                                                             %
 
 %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);
 
 %obj.volume = V_visual;
 
 obj.density = [mean_CO, mean_SC, mean_ST, mean_T1, mean_T2, mean_T3];
 %toc
 
 outputArg = 1;
 end