diff --git a/ShoulderCase/@Glenoid/Glenoid.m b/ShoulderCase/@Glenoid/Glenoid.m index 84b24c6..578541f 100644 --- a/ShoulderCase/@Glenoid/Glenoid.m +++ b/ShoulderCase/@Glenoid/Glenoid.m @@ -1,245 +1,244 @@ 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 + scapula 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; + obj.SCase = SCase; + obj.scapula = []; 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 warning(['Could not find file: ' fileName]); outputArg = 0; return; % 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: + % 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; - + 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; - + 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,:); - + % 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; + version = -atan(xProj/zProj); + inclination = -atan(yProj/zProj); 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 end end - diff --git a/ShoulderCase/@Glenoid/calcDensity.m b/ShoulderCase/@Glenoid/calcDensity.m index ab9595e..3552601 100755 --- a/ShoulderCase/@Glenoid/calcDensity.m +++ b/ShoulderCase/@Glenoid/calcDensity.m @@ -1,916 +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; +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 = (SCase.shoulder.scapula.glenoid.surface.points - origin) * R ; +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; 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); % 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 diff --git a/ShoulderCase/Acromion.m b/ShoulderCase/Acromion.m index 8458f56..c4b6868 100644 --- a/ShoulderCase/Acromion.m +++ b/ShoulderCase/Acromion.m @@ -1,213 +1,214 @@ classdef Acromion < handle %ACROMION Properties and methods associted to the acromion. % Detailed explanation goes here - + properties AI % Acromion Inddex (doi:10.2106/JBJS.D.03042) CSA % Critital Shoulder Angle (doi:10.1302/0301-620X.95B7.31028) PSA % Posterior Slope Angle (10.1016/S1058-2746(05)80036-9) PSL % Length of segment beteeen AA and AC AAA % Angulus Angle Angle (Angle between AA-origin and PA axis) AAL % Length of segment between origine and AA AAx % PA (x) position of AA AAy % IS (y) position of AA AAz % ML (z) position of AA ACx % PA (x) position of AC ACy % IS (y) position of AC ACz % ML (z) position of AC - + comment end - + properties (Hidden = true) SCase + scapula end - + methods (Access = ?Scapula) % Only Scapula is allowed to construct a Acromion function obj = Acromion(SCase) %ACROMION Construct an instance of this class % Instance of the acromion objet set all properties to zero; obj.AI = []; obj.CSA = []; obj.PSA = []; obj.PSL = []; obj.AAA = []; obj.AAL = []; obj.AAx = []; obj.AAy = []; obj.AAz = []; obj.ACx = []; obj.ACy = []; obj.ACz = []; obj.comment = ''; obj.SCase = SCase; + obj.scapula = []; end end - + methods function outputArg = morphology(obj) %MORPHOLOGY Performs morphology analysis of the acromion. % It caluculates acromion index (AI), critical shoulder angle % (CSA), and posterior slope (PS). The glenoid and scapula % surface are re-oriented (rotated and translated) in the % scapular coordinate system. For AI nd CSA, glenoid and % scapula points are projected in the scapular plane, which is % [1 0 0] after the re-orientation. % TODO: % Might be removed anatomy for concistency SCase = obj.SCase; - + %% Get local coordinate system from parent scapula - origin = SCase.shoulder.scapula.coordSys.origin; - xAxis = SCase.shoulder.scapula.coordSys.PA; - yAxis = SCase.shoulder.scapula.coordSys.IS; - zAxis = SCase.shoulder.scapula.coordSys.ML; - + origin = obj.scapula.coordSys.origin; + xAxis = obj.scapula.coordSys.PA; + yAxis = obj.scapula.coordSys.IS; + zAxis = obj.scapula.coordSys.ML; + % Rotation matrix to align to scapula coordinate system R = [xAxis' yAxis', zAxis']; - + %% Scapula, glenoid, and scapula landmarks in scapula coordinate system - + % Scapula surface alignment is done in next section, for the % caulation of AI, because of mix between auto, manual and no % segmentation. This should be corrected when there will be two % scapula object per shoulder (manual and auto). - + % Glenoid surface in scapular coordinate system - glenSurf = (SCase.shoulder.scapula.glenoid.surface.points - origin) * R; + glenSurf = (obj.scapula.glenoid.surface.points - origin) * R; % Glenoid center in scapular coordinate system - glenCenter = (SCase.shoulder.scapula.glenoid.center - origin) * R; + glenCenter = (obj.scapula.glenoid.center - origin) * R; % Scapula landmarks in scapular coordinate system - AC = (SCase.shoulder.scapula.acromioClavicular - origin) * R; % Acromio-clavicular landmark - AA = (SCase.shoulder.scapula.angulusAcromialis - origin) * R; % Angulus acromialis landmark + AC = (obj.scapula.acromioClavicular - origin) * R; % Acromio-clavicular landmark + AA = (obj.scapula.angulusAcromialis - origin) * R; % Angulus acromialis landmark obj.ACx = AC(1); obj.ACy = AC(2); obj.ACz = AC(3); obj.AAx = AA(1); obj.AAy = AA(2); obj.AAz = AA(3); - + %% Acromion Index (AI) % Adapted by AT (2018-07-18) from Valerie Mafroy Camine (2017) - + % AI = GA/GH, where: % GA: GlenoAcromial distance in scapula plane = distance from % AL to glenoid principal axis % GH: GlenoHumeral distance = 2*HHRadius, humeral head diameter (radius*2) % AL: most lateral point of the acromion - + % Get all required data, aligned in scapular plane, and % project in scapular plane. - + ScapPlaneNormal = [1 0 0]; % Normal of the scapular plane in the scapular coordinate system PlaneMean = [0 0 0]; % Origin of scapular system in scapular coordinate system - + % Project glenoid surface in scapular plane glenSurf = project2Plane(glenSurf,ScapPlaneNormal,PlaneMean,size(glenSurf,1)); - + % Project glenoid center in scapular plane glenCenter = project2Plane(glenCenter,ScapPlaneNormal,PlaneMean,size(glenCenter,1)); - + % If scapula is segmented, get AL from most lateral part of the % scapula, otherwise use acromio-clavicular landmark % Get segmentation propertiy from parent scapula - segmentedScapula = SCase.shoulder.scapula.segmentation; + segmentedScapula = obj.scapula.segmentation; segmentedScapula = strcmp(segmentedScapula,'A') || strcmp(segmentedScapula,'M'); if segmentedScapula == 1 % Rotate scapula surface to align the scapular plane with % the YZ plane, and then take the max Z (lateral) point % Transform scapula surface to scapula coordinate system - scapSurf = SCase.shoulder.scapula.surface.points; + scapSurf = obj.scapula.surface.points; scapSurf = (scapSurf - origin) * R; scapSurf = project2Plane(scapSurf,ScapPlaneNormal,PlaneMean,size(scapSurf,1)); - + % Take the most lateral point of the scapula, assuming that % this point is the most lateral point of the acromion [~,ALpositionInScapula] = max(scapSurf(:,3)); AL = scapSurf(ALpositionInScapula,:); else % No scapula points, we approximate AL with the acromioClavicular % landmark, in the scapula coordinate system AL = AC; % Acroomio-clavicalar scapula landmark AL = project2Plane(AL,ScapPlaneNormal,PlaneMean,size(AL,1)); end - + % Find glenoid principal axis with most superior and most inferior projected % glenoid points % AT: This method is not ideal. PCA would be better. It is also % used below by the CSA glenPrinAxis = [... glenSurf(glenSurf(:,2) == min(glenSurf(:,2)),:) ;... glenSurf(glenSurf(:,2) == max(glenSurf(:,2)),:)... ]; - + % Compute GA (distance from AL to glenoid principal axis) % Most inferior point of the scapula surface IG = glenPrinAxis(1, :); % Most superior point of the scapula surface SG = glenPrinAxis(2, :); GA = norm(cross(SG - IG, AL - IG))/norm(SG - IG); - + % GH (Humeral head diameter) % get humeral head radius from associated humerus GH = 2 * SCase.shoulder.humerus.radius; - + % Acromion index obj.AI = GA/GH; - + %% Critical Shoulder Angle (CSA) % Adapted by AT (2018-07-18) from Bharath Narayanan (2018) % [radCSA, degCSA] = computeCSA(GlenoidPtsinScapulaPlane, ALinScapulaPlane); - + % Vectors connecting IG to SG, and IG to AL IGSG = SG - IG; IGAL = AL - IG; CSA = vrrotvec(IGSG,IGAL); CSA = CSA(4); CSA = rad2deg(CSA); - + obj.CSA = CSA; - + %% Posterior Slope Angle and length (PSA & PSL) % By AT (2018-08-10) - + % Project AA & AC in PA-IS plane AAxy = [AA(1:2) 0]; % Juste take x and y component ACxy = [AC(1:2) 0]; % Juste take x and y component PSv = ACxy - AAxy; % Posterior slope vector ISv = [1 0 0]; % IS axis PSA = vrrotvec(PSv, ISv); PSA = PSA(4); PSA = rad2deg(PSA); - + PSL = norm(PSv); - + obj.PSA = PSA; obj.PSL = PSL; - + %% Acromial Angle Angle and length (AAA & AAL) % By AT (2018-09-13) - + % Vector between scapular origin and AA in the plane PA-IS AAv = [AA(1:2) 0]; - + % Angle between AAvect and PA axis PAv = [1 0 0]; % PS axis AAA = vrrotvec(AAv, PAv); AAA = AAA(4); AAA = rad2deg(AAA); AAA = 180 - AAA; AAL = norm(AAv); - + obj.AAA = AAA; obj.AAL = AAL; - + %% We might save an image of AI, CSA, PS - + outputArg = 1; end end end - diff --git a/ShoulderCase/Humerus.m b/ShoulderCase/Humerus.m index 89623b0..d253caf 100644 --- a/ShoulderCase/Humerus.m +++ b/ShoulderCase/Humerus.m @@ -1,258 +1,258 @@ classdef Humerus < handle %HUMERUS Properties and methods associated to the humerus % Detailed explanation goes here % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2019-06-29 % % TO DO: % Local coordinate system properties landmarks % 5 3D points center % Center of the humeral head (sphere fit on 5 points radius % Radius of the humeral head (sphere fit on 5 points jointRadius % Radius of cuvature of the articular surface (todo) SHSAngle % Scapulo-humeral subluxation angle SHSPA % Scapulo-humeral subluxation angle in the postero-anterior direction (posterior is negative, as for glenoid version) SHSIS % Scapulo-humeral subluxation angle in the infero-superior direction SHSAmpl % Scapulo-humeral subluxation (center ofset / radius) SHSOrient % Scapulo-humral subluxation orientation (0 degree is posterior) GHSAmpl % Gleno-humeral subluxation (center ofset / radius) GHSOrient % Gleno-humral subluxation orientation (0 degree is posterior) end properties (Hidden = true) SCase end methods (Access = ?Shoulder) % Only Shoulder is allowed to construct Humerus function obj = Humerus(SCase) %HUMERUS Construct an instance of this class % Detailed explanation goes here obj.landmarks = []; obj.center = []; obj.radius = []; obj.jointRadius = []; obj.SHSAngle = []; obj.SHSPA = []; obj.SHSIS = []; obj.SHSAmpl = []; obj.SHSOrient = []; obj.GHSAmpl = []; obj.GHSOrient = []; obj.SCase = SCase; end end methods function outputArg = load(obj) %LOAD Load 5 humeral head landmarks % TODO: Should be check for 'old' vs 'new' file definition SCase = obj.SCase; SCaseId4C = SCase.id4C; amiraDir = SCase.dataAmiraPath; % Try 1st with new defiunition of HHLandmarks fileName = ['newHHLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); else % Check for old version of HHLandmarks fileName = ['HHLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); warning([SCaseId4C ' is using old humeral landmarks definition']) else error('MATLAB:rmpath:DirNotFound',... 'Could not find file: %s',fileName) end end % with other methods % Check validity of data in file % Variable points belos should be renamed landmarks, for consistency try points = importedData.data; catch % Cant get landarks from file warning([SCaseId4C ' cant load humeral landmarks']) outputArg = 0; return; end % We should add a test for validity of points obj.landmarks = points; outputArg = 1; % Should report on loading result end - function outputArg = subluxation(obj) + function outputArg = subluxation(obj,scapula) % Calcul of humeral head subluxation % The function also get center and radius of the humeral head. % TODO: % Calcul would be simpler in the scapular cordinate system, as % done in Acromion object. % Radius of curvature of humeral head articular surface % Might be removed anatomy for concistency SCase = obj.SCase; %% Get data from parent shoulder object % Shoulder side side = SCase.shoulder.side; %% Get data from parent scapula object % Scapula coordinate system - xAxis = SCase.shoulder.scapula.coordSys.PA; - yAxis = SCase.shoulder.scapula.coordSys.IS; - zAxis = SCase.shoulder.scapula.coordSys.ML; + xAxis = scapula.coordSys.PA; + yAxis = scapula.coordSys.IS; + zAxis = scapula.coordSys.ML; %% Get data from parent glenoid object % Glenoid center - glenCenter = SCase.shoulder.scapula.glenoid.center; + glenCenter = scapula.glenoid.center; % Glenoid centerLine from parent glenoid object - glenCenterline = SCase.shoulder.scapula.glenoid.centerLine; + glenCenterline = scapula.glenoid.centerLine; %% Get data from humerus % humeral head landmarks (5 or ?) humHead = obj.landmarks; % Humeral head landmarks %% Anatomical measurements % Fit sphere on humeral head landmarks [humHeadCenter,humHeadRadius,HHResiduals, R2HH] = fitSphere(humHead); humHeadCenter = humHeadCenter'; % transform vertical to horizantal vector %% Scapulo-humeral subluxation (SHS) % Vector from glenoid center to humeral head center SHSvect = humHeadCenter - glenCenter; % Vector from humeral head center to scapular (z) axis, % perpendicular to z-axis --> similar to projection in xy plane % This might be rewritten more clearly SHSvectxy = SHSvect - (dot(SHSvect,zAxis)*zAxis); % SHS amplitude (ratio rather than amplitude) SHSAmpl = norm(SHSvectxy) / humHeadRadius; % SHS orientation xProj = dot(SHSvectxy, xAxis); yProj = dot(SHSvectxy, yAxis); if xProj > 0 % anterior side if yProj > 0 % ant-sup quadran SHSOrient = pi - atan(yProj/xProj); else % ant-inf quadran SHSOrient = -pi - atan(yProj/xProj); end elseif xProj < 0 % posterior side SHSOrient = - atan(yProj/xProj); else % xProj = 0 SHSOrient = pi/2 - pi * (yProj < 0); end SHSOrient = rad2deg(SHSOrient); % SHSAngle: angle between glenHumHead and zAxis SHSAngle = vrrotvec(SHSvect,zAxis); SHSAngle = SHSAngle(4); SHSAngle = rad2deg(SHSAngle); % SHSAP is angle between zxProjection and zAxis % Project glenHumHead on yAxis yProjVect = dot(SHSvectxy, yAxis)*yAxis; % Remove projection from glenHumHead --> zxProjection glenHumHeadzx = SHSvect - yProjVect; a = glenHumHeadzx; b = zAxis; rotVect4 = vrrotvec(a,b); % Unit rotation vector to align a to b rotVect3 = rotVect4(1:3); rotAngle = rotVect4(4); rotSign = sign(dot(rotVect3,yAxis)); % Sign of rotation SHSPA = rotAngle * rotSign; SHSPA = -SHSPA; % Convention that posterior is negative if strcmp(side,'L') % Change signe for left-handed scapula coordinate system SHSPA = -SHSPA; end SHSPA = rad2deg(SHSPA); % Angle in degrees % SHSIS is angle between yzProjection and zAxis % Project glenHumHead on xAxis xProjVect = dot(SHSvectxy, xAxis)*xAxis; % Remove projection from glenHumHead --> yzProjection glenHumHeadyz = SHSvect - xProjVect; a = glenHumHeadyz; b = zAxis; rotVect4 = vrrotvec(a,b); % Unit rotation vector to align a to b rotVect3 = rotVect4(1:3); rotAngle = rotVect4(4); rotSign = sign(dot(rotVect3,xAxis)); % Sign of rotation SHSIS = rotAngle * rotSign; SHSIS = SHSIS; % Convention that superior is positive if strcmp(side,'L') % Change signe for left-handed scapula coordinate system SHSIS = -SHSIS; end SHSIS = rad2deg(SHSIS); % Angle in degrees %% Angle between plane of maximum SH subluxation and CT axis [0,0,1] % Not used anymore % SHSPlane = cross(glenHumHead, zAxis); % G = vrrotvec(SHSPlane, [0 0 1]); % SHSPlaneAngle = rad2deg(G(4)); % if SHSPlaneAngle > 90 % SHSPlaneAngle = 180 - SHSPlaneAngle; % end %% Gleno-humeral subluxation (GHS) % Vector from glenoid sphere centre to glenoid centerline (perpendicular) glenCenterlineNorm = glenCenterline/norm(glenCenterline); GHS = SHSvect - dot(SHSvect, glenCenterlineNorm)*glenCenterlineNorm; % GHS amplitude GHSAmpl = norm(GHS) / humHeadRadius; % GHS orientation xProj = dot(GHS, -xAxis); yProj = dot(GHS, yAxis); if xProj ~= 0 if yProj ~= 0 GHSOrient = atan(yProj/xProj); else GHSOrient = pi * (xProj > 0); end else GHSOrient = pi/2 - pi * (yProj > 0); end GHSOrient = rad2deg(GHSOrient); %% Angle between plane of maximum GH subluxation and CT axis [0,0,1] % Not used anymore % GHSPlane = cross(glenCenterline, glenHumHead); % GHSPlaneAngle = vrrotvec(GHSPlane, [0 0 1]); % GHSPlaneAngle = rad2deg(GHSPlaneAngle(4)); % if GHSPlaneAngle > 90 % GHSPlaneAngle = 180 - GHSPlaneAngle; % end %% Set humerus properties obj.center = humHeadCenter; obj.radius = humHeadRadius; obj.SHSAngle = SHSAngle; obj.SHSPA = SHSPA; obj.SHSIS = SHSIS; obj.SHSAmpl = SHSAmpl; obj.SHSOrient = SHSOrient; obj.GHSAmpl = GHSAmpl; obj.GHSOrient = GHSOrient; % Commented ince not used % outputArgStruct = struct('HHResiduals',HHResiduals,'R2HH',R2HH); % Useful? outputArg = 1; end end end diff --git a/ShoulderCase/Scapula.m b/ShoulderCase/Scapula.m index 73d0d44..19e6a92 100644 --- a/ShoulderCase/Scapula.m +++ b/ShoulderCase/Scapula.m @@ -1,349 +1,365 @@ classdef Scapula < handle %SCAPULA Summary of this class goes here % This class defines the scapula. Landmarks are used to define its % coordinate system. It includes the glenoid object. properties angulusInferior % Landmark Angulus Inferior read from amira angulusInferiorLocal % Same as above, expressed the scapular coordinate system rather than the CT coordinate system trigonumSpinae % Landmark Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine read from amira processusCoracoideus % Landmark Processus Coracoideus (most lateral point) read from amira acromioClavicular % Landmark Acromio-clavicular joint (most lateral part on acromion)read from amira angulusAcromialis % Landmark Angulus Acromialis (most laterodorsal point) read from amira spinoGlenoidNotch % Landmark Spino-glenoid notch read from amira pillar % 5 landmarks on the pillar groove % 5 landmarks on the scapula (supraspinatus) groove coordSys % Scapular coordinate system surface % surface points and triangles of the scapula segmentation % either none 'N', manual 'M' or automatic 'A' glenoid % Glenoid class acromion % Acromion class comment % any comment end properties (Hidden = true) SCase end methods (Access = ?Shoulder) % Only Shoulder is allowed to construct a Scapula function obj = Scapula(SCase) %SCAPULA Construct an instance of this class % Constructor of the scapula object, sets all properties to % zero. obj.SCase = SCase; % All definition below should be replaced by [] obj.angulusInferior = []; obj.angulusInferiorLocal = []; obj.trigonumSpinae = []; obj.processusCoracoideus = []; obj.acromioClavicular = []; obj.angulusAcromialis = []; obj.spinoGlenoidNotch = []; obj.pillar = []; obj.groove = []; obj.surface = []; obj.segmentation = 'N'; obj.coordSys.origin = []; obj.coordSys.PA = []; % X postero-anterior axis obj.coordSys.IS = []; % Y infero-superior axis obj.coordSys.ML = []; % Z medio-lateral axis obj.glenoid = Glenoid(SCase); + obj.glenoid.scapula = obj; obj.acromion = Acromion(SCase); + obj.acromion.scapula = obj; obj.comment = ''; end end methods function outputArg = load(obj) % LOAD Load segmented surface and landmnarks % Load the segmented scapula surface (if exist) and the % scapula landmarks (6 scapula, 5 groove, 5 pillar) from % amira directory. SCase = obj.SCase; SCaseId4C = SCase.id4C; amiraDir = SCase.dataAmiraPath; matlabDir = SCase.dataMatlabPath; outputArg = 1; % Should report on loading result % Scapula (6) landmarks fileName = ['ScapulaLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); else disp(['MATLAB:rmpath:DirNotFound',... 'Could not find file: ',fileName]); outputArg = 0; return; end % Check that imported data are well formed try landmarks = importedData.data; catch outputArg = 0; return; end obj.angulusInferior = landmarks(1,:); obj.trigonumSpinae = landmarks(2,:); obj.processusCoracoideus = landmarks(3,:); obj.acromioClavicular = landmarks(4,:); obj.angulusAcromialis = landmarks(5,:); obj.spinoGlenoidNotch = landmarks(6,:); % Scapula suprapinatus groove (5) landmarks fileName = ['ScapulaGrooveLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); % Check that imported data are well formed try landmarks = importedData.data; catch outputArg = 0; return; end else disp(['MATLAB:rmpath:DirNotFound',... 'Could not find file: ',fileName]); outputArg = 0; end % Check groove landmarks validity if length(landmarks) > 5 landmarks = landmarks(1:5,:); warning(['More than 5 groove landmarks(' SCase.id ')']); elseif length(landmarks) < 5 error('Less than 5 groove landmarks'); end obj.groove = landmarks; % Scapula pillar (5) landmarks fileName = ['ScapulaPillarLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); % Check that imported data are well formed try landmarks = importedData.data; catch outputArg = 0; return; end else disp(['MATLAB:rmpath:DirNotFound',... 'Could not find file: ',fileName]); outputArg = 0; end obj.pillar = landmarks; % Import scapula surface % Try manual segmentation in amira folder fileName = ['scapula_' SCaseId4C '.stl']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 try [points,faces,~]=loadStl(fileName,1); obj.surface.points = points; obj.surface.faces = faces; obj.segmentation = 'M'; % Manual segmentation from amira catch obj.segmentation = 'E'; % Error in loading end else obj.segmentation = 'N'; % No segmentation end end function outputArg = loadAuto(obj) % Load 2 files obtained by SSM auto segmentation % scapulaSurfaceAuto{L/R}.ply % scapulaLandmarksAuto{L/R}.mat outputArg = 1; % To be set to 1 of loading is OK and 0 otherwise SCase = obj.SCase; SCaseId4C = SCase.id4C; matlabDir = SCase.dataMatlabPath; side = SCase.shoulder.side; % We might no have it yet !! + if isempty(side) + side = 'R'; + end + % Try loading auto segmentation from matlab dir fileName = ['scapulaSurfaceAuto' side '.ply']; fileName = [matlabDir '/' fileName]; + + if ~exist(fileName,'file') + side = 'L'; + fileName = ['scapulaSurfaceAuto' side '.ply']; + fileName = [matlabDir '/' fileName]; + end + + + %side = SCase.shoulder.side; + if exist(fileName,'file') == 2 try face_index_start = 1; [points,faces] = loadPly(fileName,face_index_start); % % Load ply file % ptCloud = pcread(fileName); % load the pointCloud object in the ply file % points = ptCloud.Location; % get the array of (x,y,z) points obj.surface.points = points; obj.surface.faces = faces; obj.segmentation = 'A'; % Auto segmentation from matlab catch obj.segmentation = 'E'; % Error on loading outputArg = 0; end else obj.segmentation = 'N'; % No segmentation file outputArg = 0; end % Try loading auto scapula landmarks from matlab dir filename = ['scapulaLandmarksAuto' side '.mat']; filename = [matlabDir '/' filename]; if exist(filename,'file') == 2 try % Load mat file with scapula landmarks load(filename, 'ScapulaLandmarks'); % Set loaded landmarks to scapulaAuto (obj) obj.angulusInferior = ScapulaLandmarks.angulusInferior; obj.trigonumSpinae = ScapulaLandmarks.trigonumSpinae; obj.processusCoracoideus = ScapulaLandmarks.processusCoracoideus; obj.acromioClavicular = ScapulaLandmarks.acromioClavicular; obj.angulusAcromialis = ScapulaLandmarks.angulusAcromialis; obj.spinoGlenoidNotch = ScapulaLandmarks.spinoGlenoidNotch; obj.pillar = ScapulaLandmarks.pillar; obj.groove = ScapulaLandmarks.groove; catch disp(SCaseId4C); % For debug (2 cases) outputArg = 0; end else outputArg = 0; end end function outputArg = coordSysSet(obj) % Calculate the (EPFL) scapular coordinate system % The EPFL scapular coordinate system is defined for a right scapula see % paper (10.1302/0301-620X.96B4.32641). The X axis is % antero-posterior, the Y axis is infero-superior, the Z axis % is meddio-lateral. This system is right-handed. % This orientation corresponds approximatively to the ISB % recommendadtion (doi:10.1016/j.jbiomech.2004.05.042). % For left scapula we keep the x axis as postero-anterior and % get a left-handed coordinate system. % Scapular plane is fitted on 3 points (angulusInferior, % trigonumSpinae, most laretal scapular goove landmark). % When the scapular coordinate is set, we then express the % local scapula landmarks in this system. SCase = obj.SCase; % Define local variables for landmarks aiLand = obj.angulusInferior; tsLand = obj.trigonumSpinae; pcLand = obj.processusCoracoideus; % acLand = obj.acromioClavicular; % not used here aaLand = obj.angulusAcromialis; snLand = obj.spinoGlenoidNotch; % spLand = obj.pillar; % not used here sgLand = obj.groove; % Test validity of groove if length(sgLand) ~= 5 outputArg = 0; return end % Define approximate anatomical axes (not orthogonal) postAnt = pcLand - aaLand; % ~ X infSup = snLand - aiLand; % ~ Y medLat = snLand - tsLand; % ~ Z % Find scapula side (right/left) and set to shoulder.side infSupRight = cross(medLat, postAnt); % inf-sup for right if dot(infSup, infSupRight) > 0 SCase.shoulder.side = 'R'; % Right scapula else SCase.shoulder.side = 'L'; % Left scapula end % Fit scapular plane to aiLand, tsLand and most distal scapula % groove landmarks % Determine position of most lateral landmark of the scapula groove % by calculating the distance between ts landmark and each of the 5 % landmarks of the scapula groove. dist = zeros(length(sgLand),1); for i=1:length(sgLand); dist(i) = sqrt(... (tsLand(1)-sgLand(i,1))^2 + ... (tsLand(2)-sgLand(i,2))^2 + ... (tsLand(3)-sgLand(i,3))^2 ... ); end maxDist = max(dist); % Maximum distance pos = dist == maxDist; % Location of point with max distance in the array sgLatLan = sgLand(pos,:); % Coordnates of most lateral point of the scapula groove sgLatLan = sgLatLan(1,:); % If there are more than 1, take the 1st % Fit scapular plane on aiLand, tsLand and sgLat [scapPlane, PlaneMean, ~, PlaneRMSE, ~] = ... fitPlane([aiLand; tsLand; sgLatLan]); scapPlane = scapPlane'; % Unit normal row vector % Set scapulaPlaneNormal oriented posterior-anterior (as x axis) if dot(scapPlane, postAnt) < 0 scapPlane = -scapPlane; end PAaxis = scapPlane; % PA (X) axis correspond to scapular plane normal % Scapular transverse (medio-lateral) axis % 1) Project groove points to scapular plane % 2) Fit axis to projected points % 3) Set medio-lateral direction sgProj = project2Plane(sgLand,scapPlane,sgLatLan,size(sgLand,1)); [grooveAxis, grooveMean, ~, grooveRMSE, R2Line] = fitLine(sgProj); transAxis = (grooveAxis/norm(grooveAxis))'; % Set transvere axis direction to med-lat if dot(transAxis, medLat) < 0 transAxis = -transAxis; end MLaxis = transAxis; % ML (Z) axis along scapular transverse axis % Infero-superio axis ISaxis = cross(PAaxis, MLaxis); if dot(ISaxis, infSup) < 0 % Check orientation of axis ISaxis = -ISaxis; end % !! Code below is not fully correct: grooveMean is not in % scapPlane nor in scapAxis % Set the origin of the coordinate system to the spino-gemoid notch % projected on the scapular axis transAxisMean = project2Plane(grooveMean,scapPlane,PlaneMean,1); % Project spinoGlenoidNotch to transverseAxis to set Origin origin = transAxisMean + ... dot((snLand - transAxisMean),transAxis)*transAxis; obj.coordSys.origin = origin; obj.coordSys.PA = PAaxis; obj.coordSys.IS = ISaxis; obj.coordSys.ML = MLaxis; % Express the scapula landmarks in the scapular coordinate % system R = [PAaxis' ISaxis', MLaxis']; % Rotation matrix to align to scapula coordinate system obj.angulusInferiorLocal = (obj.angulusInferior - obj.coordSys.origin) * R; % Not used % outputStuct = struct('PlaneRMSE',PlaneRMSE, 'GrooveRMSE', grooveRMSE, 'R2Line', R2Line); outputArg = 1; end end end diff --git a/measureSCase.m b/measureSCase.m index 46462c3..6c9e68e 100755 --- a/measureSCase.m +++ b/measureSCase.m @@ -1,625 +1,675 @@ function [status,message] = measureSCase(varargin) % MEASURESCASE Update the entire SCase DB with anatomical measurements. % The function can be run from (lbovenus) server with the following shell command % cd /home/shoulder/methods/matlab/database; matlab -nodisplay -nosplash -r "measureSCase;quit" % Progress can be checked through log file with following shell command % cd /home/shoulder/methods/matlab/database;tail -f log/measureSCase.log % It take ~30 min for 700 cases when run from server. % After run from server, the output file % /shoulder/data/Excel/xlsFromMatlab/anatomy.csv % must be open from Excel and saved % as xls at the same location. Excel file % /shoulder/data/Excel/ShoulderDataBase.xlsx should also be open, updated, and saved. % The script measureSCase was replacing 3 previous functions % 1) rebuildDatabaseScapulaMeasurements.m --> list of all SCases --> execute scapula_measure.m % 2) scapula_measure.m --> execute scapula_calculation and save results in Excel and SQL % 3) scapula_calculation.m --> perform anatomical analysis --> results varaiable % Input arguments can be one or several among: % {Empty,'GlenoidDensity', 'update', '[NP][1-999]', '[NP]'} % % If there is no argument the function will measure every new case. % The last 'N' or 'P' given as an argument will force the function to be % executed on all the corresponding 'N' or 'P' cases. This will overwrite % every other ShoulderCase argument '[NP][1-999]' given before. % % Examples: % measureSCase('GlenoidDensity') % is a valid instruction where every new case is measured % including its glenoid density. % measureSCase('P400','N520') % is a valid instruction where the cases 'P400' and 'N520' are % measured only if they have not already been measured yet (new % cases). % measureSCase('P400','N520','update') % is a valid instruction where the cases 'P400' and 'N520' are % measured whether they have already been measured or not. % measureSCase('GlenoidDensity','P400','N520','update','P') % is a valid instruction where all the 'P' cases are measured % including their glenoid density whether they have been % measured or not. % measureSCase('Z1590','trollingArgumentLol') % is a valid instruction which raises a warning about the % arguments and will result in the same as executing % measureSCase() i.e. every new case is calculated. % Output % File /shoulder/data/Excel/xlsFromMatlab/anatomy.csv % File /shoulder/data/Excel/xlsFromMatlab/anatomy.xls (only windows) % File /home/shoulder/data/matlab/SCaseDB.mat % File {SCaseId}.mat in each {SCaseId}/{CT}/matlab directory % logs in /log/measureSCase.log % Example: measureSCase % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2019-06-29 % TO DO: % Read first last anatomy.csv and/or SCaseDB.mat % Fill with patient and clinical data, from excel or MySQL % use argument to redo or update % Save csv within the main loop (not sure it's a good idea) % Add more message in log % Add more test to assess validity of data % Update MySQL % Rename all objects with a capital % Make it faster by not loading scapula surface if already in SCaseDB.mat % Add argument 'load' to force reload files even if they are alerady in DB % Add argument 'measure' to force re-measuring event if measurement is already in DB % Add "scapula_addmeasure.m" --> should be coded in ShoulderCase class startTime = datetime('now'); %% Open log file logFileID = openLogFile('measureSCase.log'); %% Set the data directory from the configuration file config.txt dataDir = openConfigFile('config.txt', logFileID); %% Location of the XLS ShoulderDatabase xlsDir = [dataDir '/Excel/xlsFromMatlab']; matlabDir = [dataDir '/matlab']; %% Get list of all SCase % Get list of SCases from varargin % Delault value for varargin calcGlenoidDensity = false; updateCalculations = false; specificCases = true; fprintf(logFileID, '\n\nList of SCase'); cases = {}; for argument = 1:nargin switch lower(varargin{argument}) % Test every argument in varargin. case 'glenoiddensity' calcGlenoidDensity = true; case 'update' updateCalculations = true; case regexp(lower(varargin{argument}),'[np]','match') specificCases = false; % Detect if argument is 'N' or 'P'. cases = {varargin{argument}}; % ends here. case regexp(lower(varargin{argument}),'[np][1-9][0-9]{0,2}','match') % Test for a correct argument [NP][1-999] if specificCases cases{end+1} = upper(varargin{argument}); % for element = 1:length(cases)-1 % Check if the case has already been added to the case list before if and( (upper(varargin{argument}(1))==cases{element}(1)),... % (upper(varargin{argument}(end-2:end))==cases{element}(end-2:end)) ) cases = cases(1:end-1); % Delete the last element added to cases if it is already in the list end end end otherwise warning(['\n "%s" is not a valid argument and has not been added to the'... ' list of cases.\n measureSCase arguments format must be "[npNP]",'... ' "[npNP][1-999]", "GlenoidDensity", or "update"'],varargin{argument}) end end SCaseList = listSCase(cases); %% Load current xls database (for clinical data) fprintf(logFileID, '\nLoad xls database\n\n'); addpath('XLS_MySQL_DB'); filename = [dataDir '/Excel/ShoulderDataBase.xlsx']; excelRaw = rawFromExcel(filename); % cell array excelSCaseID = excelRaw(2:end, 1); excelDiagnosis = excelRaw(2:end, 4); excelPatientHeight = excelRaw(2:end, 23); excelGlenoidWalch = excelRaw(2:end, 55); % Lines below adapted from XLS_MySQL_DB/MainExcel2XLS % Excel.diagnosisList = diagnosisListFromExcel(Excel.Raw); % Excel.treatmentList = treatmentListFromExcel(Excel.Raw); % ExcelData.Patient = patientFromExcel(excelRaw); % Structure with patient data % Excel.shoulder = shoulderFromExcel(Excel.patient, Excel.Raw); % [Excel.SCase, Excel.diagnosis, Excel.treatment, Excel.outcome, Excel.study] = SCaseFromExcel(... % Excel.shoulder, Excel.patient, Excel.diagnosisList, Excel.treatmentList, Excel.Raw); % fprintf(logFileID, ': OK'); %% Add path to ShoulderCase class addpath('ShoulderCase'); %% Instance of a ShoulderCase object if (exist('SCase','var') > 0) clear SCase; % Check for delete method end SCase = ShoulderCase; % Instanciate a ShoulderCase object SCaseDB = ShoulderCase; SCase.dataPath = dataDir; % Set dataDir for SCase %%%% definition of analysis variables % % % manual % ok = {}; +amira = {}; scapulaLoad = {}; scapulaCoord = {}; degeneration = {}; glenoidLoad = {}; glenoidMorphology = {}; glenoidDensity = {}; humerusLoad = {}; humerusSubluxation = {}; acromionMorphology = {}; % % % auto % autoOk = {}; autoScapulaLoad = {}; autoScapulaCoord = {}; %autoDegeneration = {}; autoGlenoidLoad = {}; autoGlenoidMorphology = {}; -%autoGlenoidDensity = {}; +autoGlenoidDensity = {}; %autoHumerusLoad = {}; %autoHumerusSubluxation = {}; -%autoAcromionMorphology = {}; +autoAcromionMorphology = {}; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Start loop over SCases in SCaseList nSCase = length(SCaseList); % Number of SCases for iSCaseID = 1:nSCase SCaseID = SCaseList(iSCaseID).id; SCase(iSCaseID).id = SCaseID; SCase(iSCaseID).dataPath = dataDir; % Set dataDir for SCase output = SCase(iSCaseID).path; % Set data path of this SCase percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f'); fprintf(logFileID, ['\n___________________________________\n|\n| SCaseID: ' SCaseID ' (' percentProgress '%%)']); % There are 3 parts within this SCase loop: % 1) Load & analyses manual data from amira % 2) Load & analyses auto data from matlab (Statistical Shape Model) % 3) Load clinical data from Excel database, and set them to SCase % 4) Save SCase in file SCaseID/matlab/SCase.mat % 1) Load & analyses manual data from amira if (exist([SCase(iSCaseID).dataMatlabPath '/SCase.mat'],'file') && ~updateCalculations) % New calculation of SCase.mat or update the current calculation continue end while true tic; fprintf(logFileID, '\n| Segmentation manual '); if ~output % Continue if amira dir exists in SCase fprintf(logFileID, '\n| No amira folder'); + amira{end+1} = SCaseID; break end fprintf(logFileID, '\n| Load scapula surface and landmarks'); output = SCase(iSCaseID).shoulder.scapula.load; if ~output fprintf(logFileID, ': ERROR'); scapulaLoad{end+1} = SCaseID; break end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Set coord. syst.'); output = SCase(iSCaseID).shoulder.scapula.coordSysSet; if ~output fprintf(logFileID, ': ERROR'); scapulaCoord{end+1} = SCaseID; break; end - fprintf(logFileID, '\n| Measure muscle''s degeneration'); - output = SCase(iSCaseID).shoulder.muscles.degeneration(logFileID); - - switch output % output is reversed here to differenciate different types of degeneration errors - case 1 - degeneration{1,end+1} = SCaseID; - degeneration{2,end} = output; - degeneration{3,end} = 'dicomRead() error'; - break; - case 2 - degeneration{1,end+1} = SCaseID; - degeneration{2,end} = output; - degeneration{3,end} = 'slice() error'; - break; - case 3 - degeneration{1,end+1} = SCaseID; - degeneration{2,end} = output; - degeneration{3,end} = 'slice is to small'; - break; - case 4 - degeneration{1,end+1} = SCaseID; - degeneration{2,end} = output; - degeneration{3,end} = 'crop error'; - break; - case 5 - degeneration{1,end+1} = SCaseID; - degeneration{2,end} = output; - degeneration{3,end} = 'imwrite() gave black image'; - break; - case 6 - degeneration{1,end+1} = SCaseID; - degeneration{2,end} = output; - degeneration{3,end} = 'automatic segmentation error'; - break; - end - fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Load glenoid surface'); output = SCase(iSCaseID).shoulder.scapula.glenoid.load; if ~output fprintf(logFileID, '\n| glenoid surface loading error'); glenoidLoad{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Measure glenoid anatomy'); output = SCase(iSCaseID).shoulder.scapula.glenoid.morphology; if ~output fprintf(logFileID, ': ERROR'); glenoidMorphology{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); if calcGlenoidDensity fprintf(logFileID, '\n| Measure glenoid density'); output = SCase(iSCaseID).shoulder.scapula.glenoid.calcDensity; if ~output fprintf(logFileID, ': Density calculus error. Could not read dicom'); glenoidDensity{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); end - + fprintf(logFileID, '\n| Load humerus data'); output = SCase(iSCaseID).shoulder.humerus.load; if ~output fprintf(logFileID, ': ERROR'); humerusLoad{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Measure humerus subluxation'); - output = SCase(iSCaseID).shoulder.humerus.subluxation; + output = SCase(iSCaseID).shoulder.humerus.subluxation(SCase(iSCaseID).shoulder.scapula); if ~output fprintf(logFileID, ': ERROR'); humerusSubluxation{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Measure acromion anatomy'); % Should be run after SCase.shoulder.humerus (for AI) output = SCase(iSCaseID).shoulder.scapula.acromion.morphology; if ~output fprintf(logFileID, ': ERROR'); acromionMorphology{end+1} = SCaseID; break; end + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure muscle''s degeneration'); + output = SCase(iSCaseID).shoulder.muscles.degeneration(logFileID); + + switch output % output is reversed here to differenciate different types of degeneration errors + case 1 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'dicomRead() error'; + break; + case 2 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'slice() error'; + break; + case 3 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'slice is to small'; + break; + case 4 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'crop error'; + break; + case 5 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'imwrite() gave black image'; + break; + case 6 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'automatic segmentation error'; + break; + end + fprintf(logFileID, ': OK'); ok{end+1} = SCaseID; % 2) Set clinical data from Excel database to SCase fprintf(logFileID, '\n| Set clinical data from Excel'); % Get idx of excelRaw for SCaseID % Use cell2struct when dots are removed from excel headers idx = find(strcmp(excelRaw, SCaseID)); % Index of SCaseID in excelRaw SCase(iSCaseID).diagnosis = excelRaw{idx,4}; SCase(iSCaseID).treatment = excelRaw{idx,9}; SCase(iSCaseID).patient.gender = excelRaw{idx,19}; SCase(iSCaseID).patient.age = excelRaw{idx,21}; SCase(iSCaseID).patient.height = excelRaw{idx,23}; SCase(iSCaseID).patient.weight = excelRaw{idx,24}; SCase(iSCaseID).patient.BMI = excelRaw{idx,25}; SCase(iSCaseID).shoulder.scapula.glenoid.walch = excelRaw{idx,55}; + break; + % Patient data --> TO DO % find patient id (in SQL) corresponding to iSCaseID % % % idxSCase = find(contains(ExcelData, SCaseID)); % SCase(iSCaseID).patient.id = ExcelData(idxSCase).patient.patient_id; % SQL patient id % SCase(iSCaseID).patient.idMed = ExcelData(idxSCase).patient.IPP; % IPP to be replace by coded (anonymized) IPP from SecuTrial % SCase(iSCaseID).patient.gender = ExcelData(idxSCase).patient.gender; % SCase(iSCaseID).patient.age = []; % Age (years) at treatment, or preop CT (we might also add birthdate wi day et at 15 of the month) % SCase(iSCaseID).patient.ethnicity = ExcelData(idxSCase).patient.IPP; % SCase(iSCaseID).patient.weight = ExcelData(idxSCase).patient.weight; % SCase(iSCaseID).patient.height = ExcelData(idxSCase).patient.height; % SCase(iSCaseID).patient.BMI = ExcelData(idxSCase).pateint.BMI; % SCase(iSCaseID).patient.comment = ExcelData(idxSCase).patient.comment; - % 3) Save SCase in file SCaseID/matlab/SCase.mat - fprintf(logFileID, '\n| Save SCase.mat'); - output = SCase(iSCaseID).saveMatlab; - if output - fprintf(logFileID, ': OK'); - else - fprintf(logFileID, ': ERROR'); - end - break; end fprintf(logFileID, '\n| Elapsed time (s): %s\n|', num2str(toc)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if (exist([SCase(iSCaseID).dataMatlabPath '/scapulaLandmarksAutoR.mat'],'file') && ~updateCalculations) continue end while true tic; - fprintf(logFileID, '\n| Segmentation auto '); fprintf(logFileID, '\n| Load scapula surface and landmarks'); output = SCase(iSCaseID).shoulder.scapulaAuto.loadAuto; % Load auto data from matlab directory if ~output fprintf(logFileID, ': ERROR'); autoScapulaLoad{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Set coord. syst.'); output = SCase(iSCaseID).shoulder.scapulaAuto.coordSysSet; if ~output fprintf(logFileID, ': ERROR'); autoScapulaCoord{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Load glenoid surface'); output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.loadAuto; if ~output fprintf(logFileID, ': ERROR'); autoGlenoidLoad{end+1} = SCaseID; break; end fprintf(logFileID, ': OK'); fprintf(logFileID, '\n| Measure glenoid anatomy'); output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.morphology; if ~output fprintf(logFileID, ': ERROR'); autoGlenoidMorphology{end+1} = SCaseID; break; end + fprintf(logFileID, ': OK'); + if calcGlenoidDensity + fprintf(logFileID, '\n| Measure glenoid density'); + output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.calcDensity; + + if ~output + fprintf(logFileID, ': Density calculus error. Could not read dicom'); + autoGlenoidDensity{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + end + + fprintf(logFileID, '\n| Load humerus data'); + output = SCase(iSCaseID).shoulder.humerus.load; + + if ~output + fprintf(logFileID, ': ERROR'); + if humerusLoad{end} ~= SCaseID; + humerusLoad{end+1} = SCaseID; + end + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure humerus subluxation'); + output = SCase(iSCaseID).shoulder.humerus.subluxation(SCase(iSCaseID).shoulder.scapulaAuto); + + if ~output + fprintf(logFileID, ': ERROR'); + humerusSubluxation{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure acromion anatomy'); + % Should be run after SCase.shoulder.humerus (for AI) + output = SCase(iSCaseID).shoulder.scapulaAuto.acromion.morphology; + + if ~output + fprintf(logFileID, ': ERROR'); + autoAcromionMorphology{end+1} = SCaseID; + break; + end + fprintf(logFileID, ': OK'); autoOk{end+1} = SCaseID; + fprintf(logFileID, '\n| Save SCase.mat'); + output = SCase(iSCaseID).saveMatlab; + + if ~output + fprintf(logFileID, ': ERROR'); + break; + end + + fprintf(logFileID, ': OK'); break; end fprintf(logFileID, '\n| Elapsed time (s): %s\n|___________________________________\n\n', num2str(toc)); % 4) Load & analyses auto data from matlab (Statistical Shape Model) % Load scapula surface and landmarks save('measureSCase_analysis.mat',... - 'ok','scapulaLoad','scapulaCoord',... + 'ok','amira','scapulaLoad','scapulaCoord',... 'degeneration','glenoidLoad','glenoidMorphology','glenoidDensity',... 'humerusLoad','humerusSubluxation','acromionMorphology',... 'autoOk','autoScapulaLoad','autoScapulaCoord','autoGlenoidLoad',... - 'autoGlenoidMorphology') + 'autoGlenoidMorphology','autoGlenoidDensity','autoAcromionMorphology') end % End of loop on SCaseList %% Write csv file % This might be a function (SCase.csvSave()). The input would be a filename and a structure % data % Replace header and data by structure. Currently not working %{ txtFilename = [xlsDir, '/anatomy.txt']; % Name of the txt file DataStruc = struc(... 'SCase_id', SCase.id, ... 'shoulder_side', SCase.shoulder.side,... 'glenoid_radius', SCase.shoulder.scapula.glenoid.radius,... 'glenoid_sphereRMSE', SCase.shoulder.scapula.glenoid.sphereRMSE,... 'glenoid_depth', SCase.shoulder.scapula.glenoid.depth,... 'glenoid_width', SCase.shoulder.scapula.glenoid.width,... 'glenoid_height', SCase.shoulder.scapula.glenoid.height,... 'glenoid_centerPA', SCase.shoulder.scapula.glenoid.centerLocal(1),... 'glenoid_centerIS', SCase.shoulder.scapula.glenoid.centerLocal(2),... 'glenoid_centerML', SCase.shoulder.scapula.glenoid.centerLocal(3),... 'glenoid_versionAmpl', SCase.shoulder.scapula.glenoid.versionAmpl,... 'glenoid_versionOrient', SCase.shoulder.scapula.glenoid.versionOrient,... 'glenoid_version', SCase.shoulder.scapula.glenoid.version,... 'glenoid_inclination', SCase.shoulder.scapula.glenoid.inclination,... 'humerus_jointRadius', SCase.shoulder.humerus.jointRadius,... 'humerus_headRadius', SCase.shoulder.humerus.radius,... 'humerus_GHSAmpl', SCase.shoulder.humerus.GHSAmpl,... 'humerus_GHSOrient', SCase.shoulder.humerus.GHSOrient,... 'humerus_SHSAmpl', SCase.shoulder.humerus.SHSAmpl,... 'humerus_SHSOrient', SCase.shoulder.humerus.SHSOrient,... 'humerus_SHSAngle', SCase.shoulder.humerus.SHSAngle,... 'humerus_SHSPA', SCase.shoulder.humerus.SHSPA,... 'humerus_SHSIS', SCase.shoulder.humerus.SHSIS,... 'acromion_AI', SCase.shoulder.scapula.acromion.AI,... 'acromion_CSA', SCase.shoulder.scapula.acromion.CSA,... 'acromion_PS', SCase.shoulder.scapula.acromion.PS... ); DataTable = strct2table(Data); writetable(DataTable,filename); %} % Header of the csv file dataHeader = [... 'SCase_id,' ... 'shoulder_side,' ... 'glenoid_radius,' ... 'glenoid_sphereRMSE,' ... 'glenoid_depth,' ... 'glenoid_width,' ... 'glenoid_height,' ... 'glenoid_centerPA,' ... 'glenoid_centerIS,' ... 'glenoid_centerML,' ... 'glenoid_versionAmpl,' ... 'glenoid_versionOrient,' ... 'glenoid_version,' ... 'glenoid_inclination,' ... 'humerus_jointRadius,' ... 'humerus_headRadius,' ... 'humerus_GHSAmpl,' ... 'humerus_GHSOrient,' ... 'humerus_SHSAmpl,' ... 'humerus_SHSOrient,' ... 'humerus_SHSAngle,' ... 'humerus_SHSPA,' ... 'humerus_SHSIS,' ... 'acromion_AI,' ... 'acromion_CSA,' ... 'acromion_PSA,'... 'acromion_AAA\n'... ]; updateList = listSCase; % The list of SCase to use to construct anatomy.csv and SCaseDB.mat fprintf(logFileID, '\n\nSave anatomy.csv file'); fid = fopen([xlsDir, '/anatomy.csv'],'w'); %Last csv is deleted and a new one, updated is being created fprintf(fid,dataHeader); fclose(fid); % The following lines re-create the SCaseDB.mat and contruct the anatomy.csv file S = ShoulderCase; for Case=1:length(updateList) S(Case).id = updateList(Case).id; % S(Case).dataPath = dataDir; % S(Case).path; % Compute the SCase.dataMatlabPath loaded = load([S(Case).dataMatlabPath '/SCase.mat']); % Load SCase.mat SCaseDB(Case) = loaded.SCase; SCaseDB(Case).appendToCSV('anatomy.csv'); end fprintf(logFileID, ': OK'); %% Save the entire SCaseDB array as a matlab file fprintf(logFileID, '\n\nSave SCase database'); filename = 'SCaseDB'; filename = [matlabDir '/' filename '.mat']; try save(filename, 'SCaseDB'); fprintf(logFileID, ': OK'); catch fprintf(logFileID, ': Error'); end %fprintf(logFileID, [': ' csvFilename]); %% Write xls file (Windows only) % [~,message] = csv2xlsSCase(csvFilename); % fprintf(logFileID, message); % If run from non-windows system, only csv will be produced. The xls can be % produced by opening the csv from Excel and save as xls. % Could be a function with 1 input (cvsFilenames %{ xlsFilename = strrep(csvFilename, '.csv', '.xls'); % Replace .csv by .xls if ispc % Check if run from windows! % Read cvs file as table and write the table as xls cvsTable = readtable(csvFilename); % Get table from csv (only way to read non-numeric data) cvsCell = table2cell(cvsTable); % Tranform table cell array dataHeader = cvsTable.Properties.VariableNames; % Get header cvsCell = [dataHeader; cvsCell]; % Add header sheet = 'anatomy'; % Sheet name of the xls file [status,message] = xlswrite(xlsFilename, cvsCell, sheet); % Save xls if status fprintf(logFileId, ['\nSaved ' xlsFilename]); else fprintf(logFileId, ['\nCould not save ' xlsFilename]); fprintf(logFileId, ['\n' message.identifier]); fprintf(logFileId, ['\n' message.message]); end else fprintf(logFileId, ['\n' xlsFilename ' not saved. Open and save as xls from Excel']); end %} %% Close log file fprintf(logFileID, '\n\nSCase measured: %s', num2str(length(SCaseList))); % Number of SCases measured elapsedTime = char(datetime('now')-startTime); elapsedTime = [elapsedTime(1:2) ' hours ' elapsedTime(4:5) ' minutes ' elapsedTime(7:8) ' seconds']; fprintf(logFileID, ['\nTotal Elapsed time: ' elapsedTime ]); fprintf(logFileID, '\n'); fclose(logFileID); % Close log file %% Output of the SCase if required by argument % SCase.output % inputArg = abaqus/density/references/ % update mySQL % SCase.sqlSave status = 1; message = 'OK'; end