Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (master) Browse

measureDensity.m
5a906fe7379dmaster

measureDensity.m
View Options

			function outputArg = measureDensity(obj, varargin)
			%% 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

			Timer.start();
			Logger.log("Glenoid density: ");


			SCase = obj.scapula.shoulder.SCase;
			glenSurf = obj.surface.points;

			% Get local coordinate system from parent scapula object
			origin = obj.scapula.coordSys.origin;
			xAxis = obj.scapula.coordSys.PA;
			yAxis = obj.scapula.coordSys.IS;
			zAxis = obj.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 = (obj.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;
			Logger.logn("failed after %s", Timer.stop());
			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(round(size_list(1)/2)).name] ;

			try
			dicom_information = dicominfo(dicominfopath);
			catch
			outputArg = false;
			Logger.logn("failed after %s", Timer.stop());
			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 (j,i,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
			%(j,i,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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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(2)) && (j<=size_V(1)) && (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


			% TRY/CATCH added by AT to avoid error with P255
			try
			V_visual = V+1000*V_cyl; %can be used in volumeViewer for rendering
			catch
			outputArg = 0;
			Logger.logn("failed after %s", Timer.stop());
			return
			end

			%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

			Logger.logn("OK %s", Timer.stop());
			outputArg = 1;
			end

R8218/ShoulderCase/@Glenoid/measureDensity.m5a906fe7379dmaster

measureDensity.mView Options

R8218/ShoulderCase/@Glenoid/measureDensity.m
5a906fe7379dmaster

measureDensity.m
View Options