diff --git a/ShoulderCase/@Glenoid/Glenoid.m b/ShoulderCase/@Glenoid/Glenoid.m index d4c5a74..abc30b2 100644 --- a/ShoulderCase/@Glenoid/Glenoid.m +++ b/ShoulderCase/@Glenoid/Glenoid.m @@ -1,246 +1,242 @@ classdef Glenoid < handle %GLENOID Summary of this class goes here % Detailed explanation goes here % Need to load surface points from Amira (or auto-segment) % Need to set anatomical calculation, and use the scapular coord.Syst % TODO % loadMatlab() properties surface % points of the glenoid surface load glenoid surface center % Glenoid center in CT coordinate system centerLocal % Glenoid center in scapula coordinate system radius sphereRMSE centerLine depth width height versionAmpl versionOrient version inclination density walch comment end properties (Hidden = true) SCase end methods (Access = ?Scapula) % Only Scapula is allowed to construct a Glenoid function obj = Glenoid(SCase) %GLENOID Construct an instance of this class % Detailed explanation goes here obj.surface = []; obj.center = []; obj.centerLocal.x = []; obj.centerLocal.y = []; obj.centerLocal.z = []; obj.radius = []; obj.sphereRMSE = []; obj.depth = []; obj.width = []; obj.height = []; obj.versionAmpl = []; obj.versionOrient = []; obj.version = []; obj.inclination = []; obj.density = []; obj.walch = []; obj.comment = []; obj.SCase = SCase; end end methods function outputArg = load(obj) % LOAD Summary of this method goes here % Load the points of the glenoid surface from the amira % directoty. SCase = obj.SCase; SCaseId4C = SCase.id4C; amiraDir = SCase.dataAmiraPath; % Import glenoid surface points fileName = ['ExtractedSurface' SCaseId4C '.stl']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 [points,faces,~]=loadStl(fileName,1); obj.surface.points = points; obj.surface.faces = faces; else error(['Could not find file: ' fileName]) end outputArg = 1; % Should report on loading result end function outputArg = loadAuto(obj) % LOADAUTO Load genoid surface from auto segmentation % Load the points of the glenoid surface from matlab dir SCase = obj.SCase; SCaseId4C = SCase.id4C; matlabDir = SCase.dataMatlabPath; side = SCase.shoulder.side; % Import glenoid surface points fileName = ['glenoidSurfaceAuto' side '.ply']; fileName = [matlabDir '/' fileName]; if exist(fileName,'file') == 2 [points,faces]=loadPly(fileName,1); obj.surface.points = points; obj.surface.faces = faces; else outputArg = 0; return % error(['Could not find file: ' fileName]) end outputArg = 1; % Should report on loading result end function outputArg = morphology(obj) % Calculate morphological properties of the glenoid from its surface: % center, radius, depth, width, height, % versionAmpl, versionOrient, version, inclination 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,:); % Fit sphere on glenoid surface [sphereCenter,sphereRadius,sphereResiduals, R2Glenoid] = fitSphere(glenSurf); sphereCenter = sphereCenter'; % transform a vertical to horizotal vector glenRadius = sphereRadius; % Root Mean Square Error sphereRMSE = norm(sphereResiduals) / sqrt(length(sphereResiduals)); % Glenoid centerline centerLine = sphereCenter - glenCenter; % Depth of the glenoid surface % Project the points of the surface on the centerline glenSurfProj = zeros(size(glenSurf)); for i = 1:size(glenSurf,1) glenSurfProj(i,:) = glenCenter + (dot(centerLine,(glenSurf(i,:) - glenCenter)) / ... dot(centerLine,centerLine)) * centerLine; dist(i) = norm(glenCenter - glenSurfProj(i,:)); end depth = max(dist); % Glenoid width and height % This part should be rewritten width = ''; height = ''; % Glenoid version 3D, defined by the amplitude and orientation of the angle % between the glenoid centerline and the scapular transverse axis (z axis). % The orientation is the angle between the x axis and the glenoid % centerline projected in the xy plane. Zero orientation correspond to -x % axis (posterior side), 90 to superior, 180 to anterior orientaion, and % -90 to inferior. Glenoid version is also reported as version2D (>0 for % retroversion) and inclination (>0 for superior tilt). versionAmplRot = vrrotvec(zAxis,centerLine); % angle between centerLine and scapular axis versionAmpl = rad2deg(versionAmplRot(4)); xProj = dot(centerLine, xAxis); yProj = dot(centerLine, yAxis); zProj = dot(centerLine, zAxis); if xProj > 0 % anterior side if yProj > 0 % ant-sup quadran versionOrient = pi - atan(yProj/xProj); else % ant-inf quadran versionOrient = -pi - atan(yProj/xProj); end elseif xProj < 0 % posterior side versionOrient = - atan(yProj/xProj); else % xProj = 0 versionOrient = pi/2 - pi * (yProj < 0); end versionOrient = rad2deg(versionOrient); if zProj > 0 version = atan(xProj/zProj); inclination = atan(yProj/zProj); else % not realistic situation version = NaN; inclination = NaN; end version = rad2deg(version); inclination = rad2deg(inclination); % Angle between plane of maximum wear and CT axis [0,0,1] maxWearPlane = cross(centerLine, zAxis); wearPlaneAngle = vrrotvec(maxWearPlane, [0 0 1]); wearPlaneAngle = rad2deg(wearPlaneAngle(4)); if wearPlaneAngle > 90 wearPlaneAngle = 180 - wearPlaneAngle; end % Glenoid center in scapula coordinate system CT2scap = [xAxis' yAxis' zAxis']; glenCenterLocal = (glenCenter - origin) * CT2scap; if length(glenCenterLocal) ~= 3 warning(['Variable glenCenterLocal not defined (' SCase.id ')']); outputArg = 0; return end % Set object properties obj.center = glenCenter; obj.centerLocal.x = glenCenterLocal(1); obj.centerLocal.y = glenCenterLocal(2); obj.centerLocal.z = glenCenterLocal(3); obj.radius = glenRadius; obj.sphereRMSE = sphereRMSE; obj.depth = depth; obj.centerLine = centerLine; obj.width = width; obj.height = height; obj.versionAmpl = versionAmpl; obj.versionOrient = versionOrient; obj.version = version; obj.inclination = inclination; outputArg = 1; end - - function outputArg = calcDensity(obj) - outputArg = 1; - end end end diff --git a/ShoulderCase/@Glenoid/calcDensity.m b/ShoulderCase/@Glenoid/calcDensity.m index d4c5a74..8487779 100644 --- a/ShoulderCase/@Glenoid/calcDensity.m +++ b/ShoulderCase/@Glenoid/calcDensity.m @@ -1,246 +1,903 @@ -classdef Glenoid < handle - %GLENOID Summary of this class goes here - % Detailed explanation goes here - % Need to load surface points from Amira (or auto-segment) - % Need to set anatomical calculation, and use the scapular coord.Syst - - % TODO - % loadMatlab() - - properties - surface % points of the glenoid surface load glenoid surface - center % Glenoid center in CT coordinate system - centerLocal % Glenoid center in scapula coordinate system - radius - sphereRMSE - centerLine - depth - width - height - versionAmpl - versionOrient - version - inclination - density - walch - comment - end - - properties (Hidden = true) - SCase - end - - methods (Access = ?Scapula) % Only Scapula is allowed to construct a Glenoid - function obj = Glenoid(SCase) - %GLENOID Construct an instance of this class - % Detailed explanation goes here - obj.surface = []; - obj.center = []; - obj.centerLocal.x = []; - obj.centerLocal.y = []; - obj.centerLocal.z = []; - obj.radius = []; - obj.sphereRMSE = []; - obj.depth = []; - obj.width = []; - obj.height = []; - obj.versionAmpl = []; - obj.versionOrient = []; - obj.version = []; - obj.inclination = []; - obj.density = []; - obj.walch = []; - obj.comment = []; - - obj.SCase = SCase; +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 + + +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); + + +%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 - methods - function outputArg = load(obj) - % LOAD Summary of this method goes here - % Load the points of the glenoid surface from the amira - % directoty. - - SCase = obj.SCase; - SCaseId4C = SCase.id4C; - - amiraDir = SCase.dataAmiraPath; - - % Import glenoid surface points - fileName = ['ExtractedSurface' SCaseId4C '.stl']; - fileName = [amiraDir '/' fileName]; - if exist(fileName,'file') == 2 - [points,faces,~]=loadStl(fileName,1); - obj.surface.points = points; - obj.surface.faces = faces; - else - error(['Could not find file: ' fileName]) +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 - outputArg = 1; % Should report on loading result - end - - function outputArg = loadAuto(obj) - % LOADAUTO Load genoid surface from auto segmentation - % Load the points of the glenoid surface from matlab dir - - SCase = obj.SCase; - SCaseId4C = SCase.id4C; - - matlabDir = SCase.dataMatlabPath; - side = SCase.shoulder.side; - - % Import glenoid surface points - fileName = ['glenoidSurfaceAuto' side '.ply']; - fileName = [matlabDir '/' fileName]; - if exist(fileName,'file') == 2 - [points,faces]=loadPly(fileName,1); - obj.surface.points = points; - obj.surface.faces = faces; +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 - outputArg = 0; - return - % error(['Could not find file: ' fileName]) + continue end - - outputArg = 1; % Should report on loading result - end - - function outputArg = morphology(obj) - % Calculate morphological properties of the glenoid from its surface: - % center, radius, depth, width, height, - % versionAmpl, versionOrient, version, inclination - - 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 + 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 - glenCenterNode = find(dist == min(dist)); - glenCenter = glenSurf(glenCenterNode,:); - - % Fit sphere on glenoid surface - [sphereCenter,sphereRadius,sphereResiduals, R2Glenoid] = fitSphere(glenSurf); - sphereCenter = sphereCenter'; % transform a vertical to horizotal vector - glenRadius = sphereRadius; - % Root Mean Square Error - sphereRMSE = norm(sphereResiduals) / sqrt(length(sphereResiduals)); - - % Glenoid centerline - centerLine = sphereCenter - glenCenter; - - % Depth of the glenoid surface - % Project the points of the surface on the centerline - glenSurfProj = zeros(size(glenSurf)); - for i = 1:size(glenSurf,1) - glenSurfProj(i,:) = glenCenter + (dot(centerLine,(glenSurf(i,:) - glenCenter)) / ... - dot(centerLine,centerLine)) * centerLine; - dist(i) = norm(glenCenter - glenSurfProj(i,:)); + 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 - depth = max(dist); - - - % Glenoid width and height - % This part should be rewritten - - width = ''; - height = ''; - - - % Glenoid version 3D, defined by the amplitude and orientation of the angle - % between the glenoid centerline and the scapular transverse axis (z axis). - % The orientation is the angle between the x axis and the glenoid - % centerline projected in the xy plane. Zero orientation correspond to -x - % axis (posterior side), 90 to superior, 180 to anterior orientaion, and - % -90 to inferior. Glenoid version is also reported as version2D (>0 for - % retroversion) and inclination (>0 for superior tilt). - - versionAmplRot = vrrotvec(zAxis,centerLine); % angle between centerLine and scapular axis - versionAmpl = rad2deg(versionAmplRot(4)); - - xProj = dot(centerLine, xAxis); - yProj = dot(centerLine, yAxis); - zProj = dot(centerLine, zAxis); - - if xProj > 0 % anterior side - if yProj > 0 % ant-sup quadran - versionOrient = pi - atan(yProj/xProj); - else % ant-inf quadran - versionOrient = -pi - atan(yProj/xProj); - end - elseif xProj < 0 % posterior side - versionOrient = - atan(yProj/xProj); - else % xProj = 0 - versionOrient = pi/2 - pi * (yProj < 0); + 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 - versionOrient = rad2deg(versionOrient); - if zProj > 0 - version = atan(xProj/zProj); - inclination = atan(yProj/zProj); - else % not realistic situation - version = NaN; - inclination = NaN; + 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 - version = rad2deg(version); - inclination = rad2deg(inclination); - - % Angle between plane of maximum wear and CT axis [0,0,1] - maxWearPlane = cross(centerLine, zAxis); - wearPlaneAngle = vrrotvec(maxWearPlane, [0 0 1]); - wearPlaneAngle = rad2deg(wearPlaneAngle(4)); - if wearPlaneAngle > 90 - wearPlaneAngle = 180 - wearPlaneAngle; + 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 - - % Glenoid center in scapula coordinate system - CT2scap = [xAxis' yAxis' zAxis']; - glenCenterLocal = (glenCenter - origin) * CT2scap; - if length(glenCenterLocal) ~= 3 - warning(['Variable glenCenterLocal not defined (' SCase.id ')']); - outputArg = 0; - return + 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 - - % Set object properties - obj.center = glenCenter; - obj.centerLocal.x = glenCenterLocal(1); - obj.centerLocal.y = glenCenterLocal(2); - obj.centerLocal.z = glenCenterLocal(3); - obj.radius = glenRadius; - obj.sphereRMSE = sphereRMSE; - obj.depth = depth; - obj.centerLine = centerLine; - obj.width = width; - obj.height = height; - obj.versionAmpl = versionAmpl; - obj.versionOrient = versionOrient; - obj.version = version; - obj.inclination = inclination; - - outputArg = 1; + 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 - function outputArg = calcDensity(obj) - outputArg = 1; - end - end +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