Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- clear
- clc
- % Allow for user input to ask for which bone test case is desired
- select_testCase = input(['Select bone test case:\n\n - Homogenous Young Bone [1]'...
- '\n - Homogenous Old Bone [2]'...
- '\n - Heterogenous Young Bone [3]'...
- '\n\n >> Enter Input Value: '],'s');
- clc
- % Allow for user input to ask for implant location if desired
- select_locationImplant = input(['Select implant location (CASE SENSITIVE!):\n\n - Point A (28.08mm,38.22mm) [A]'...
- '\n - Point B (26.52mm,28.86mm) [B]'...
- '\n - Point C (32.76mm,31.98mm) [C]'...
- '\n - Point D (39.00mm,25.74mm) [D]'...
- '\n - No Implant! [N]'...
- '\n\n >> Enter Input Value: ',],'s');
- clc
- % Accquire GS values from provided 'Femur.png'
- A = imread('Femur','png');
- A = rgb2gray(A); % Convert the image to grayscale
- A = double(A); % Convert the uint8 (8-bit) to a double
- % This makes operations with A have double precision
- res = 0.78e-3; % resolution in meters per pixel
- % Assigning GS for homogenous cases (young & old/A & B)
- if strcmp(select_testCase,'1') == 1 || strcmp(select_testCase,'2') == 1
- for i=1:86
- for j=1:86
- if A(i,j) < 255/2 % Only consider BMD > 1g/cc ... BMD = 2/255 GS
- A(i,j) = 0; % if below GS threshoold, BMD = 0 (as if nothing is there!)
- end
- end
- end
- end
- % Creating Young's modulus matrix correlating each pixel of A to a modulus
- E = zeros(86); %preallocating modulus matrix, E
- switch select_testCase
- case '1' % homogenous YOUNG bone
- for i=1:86
- for j=1:86
- if A(i,j) ~= 0 %if there is a GS/BMD value then assign E
- E(i,j) = 17E9; % homogenous YOUNG bone E = 17GPa
- end
- end
- end
- case '2' % homogenous OLD bone
- for i=1:86
- for j=1:86
- if A(i,j) ~= 0 %if there is a GS/BMD value then assign E
- E(i,j) = 13E9; % homogenous OLD bone E = 13GPa
- end
- end
- end
- case '3' % hetergenous young bone
- % E = a * (BMD^b) <- power relationship
- b = (log(1.107/17)) / (log(0.26/2));
- a = 17 / (2^b);
- for i=1:86
- for j=1:86
- E(i,j) = 1E9 * a * ((2*A(i,j)/255)^b); % 1E9 to turn GPa to Pa
- end
- end
- end
- %Phase III only runs IF there is any implant location selected
- if strcmp(select_locationImplant, 'N') ~= 1
- % Accquire GS values from 'Implant.png'
- B = imread('Implant','png');
- B = rgb2gray(B); % Convert the image to grayscale
- B = double(B); % Convert the uint8 (8-bit) to a double
- % This makes operations with A have double precision
- % IF a pre-defined implant location is selected, it will assign the
- % specified location a pixel coordinate for the implant center
- switch select_locationImplant
- case 'A'
- implantCenter_pixel = [50, 37];
- case 'B'
- implantCenter_pixel = [37, 35];
- case 'C'
- implantCenter_pixel = [41, 43];
- case 'D'
- implantCenter_pixel = [34, 51];
- end
- % Assigning modulus of steel implant 200GPa to implant in the femur
- for i=1:15
- for j=1:15
- if B(i,j) ~= 0 % for areas in implant with >0 GS value
- E(implantCenter_pixel(1) - 8 + i, implantCenter_pixel(2) - 8 + j) = 200E9;
- end
- end
- end
- end
- % Find BMD weighted centroid
- sumy = 0; % preallocating
- sumz = 0;
- total_sum = 0;
- position = @(i,j) [(j-1)*res+res/2, (i-1)*res+res/2]; % location of each pixel in m (Z & Y position)
- for i = 1:86
- for j = 1:86
- posEqn = position(i,j); % vector to determine the CENTER of the pixel given pixel resolution 0.78e-3m/pixel
- sumz = sumz + posEqn(1) * E(i,j) * res^2;
- sumy = sumy + posEqn(2) * E(i,j) * res^2;
- total_sum = total_sum + E(i,j) * res^2;
- end
- end
- % Centroid coodinates
- y_hat = sumy/total_sum; % in m
- z_hat = sumz/total_sum; % in m
- centroid_pixelPosition = [z_hat/res, y_hat/res]; %turning m back into pixels
- % Given moment values in the problem statement
- M_y = 46; %in N*m
- M_z = 28; %in N*m
- % Determine moment of inertia (MoI)
- I_zz_star = 0; %preallocation
- I_yy_star = 0;
- I_yz_star = 0;
- for i=1:86
- for j=1:86
- posEqn = position(i,j); % position of pixel CENTER in m
- I_zz_star = I_zz_star + E(i,j) * (posEqn(2) - y_hat)^2 * res^2 + (res^4 / 12); % because parallel axis theorem
- I_yy_star = I_yy_star + E(i,j) * (posEqn(1) - z_hat)^2 * res^2 + (res^4 / 12);
- I_yz_star = I_yz_star + E(i,j) * (posEqn(2) - y_hat) * (posEqn(1) - z_hat) * res^2;
- end
- end
- % Determine stress distrubtion (sigmaMatrix) using previous MoI
- sigmaMatrix = zeros(86); %preallocation
- for i=1:86
- for j=1:86
- %if E(i,j) ~= 200E9 % ignores stresses from the IMPLANT
- posEqn = position(i,j); % position of pixel CENTER in m
- sigmaMatrix(i,j) = E(i,j)*((M_y*I_zz_star + M_z*I_yz_star)*(posEqn(1)-z_hat)-(M_y*I_yz_star + M_z*I_yy_star)*(posEqn(2)-y_hat))/((I_yy_star*I_zz_star)-(I_yz_star^2));
- %end
- end
- end
- max(max(sigmaMatrix))
- [row,col] = find(sigmaMatrix == max(sigmaMatrix(:)))
- z = (col-1)*.78+(.78/2)
- y = (row-1)*.78+(.78/2)
- y_hat
- z_hat
- save('sigmaMatrix')
- contourf(sigmaMatrix)
- c = colorbar;
- %caxis([-1E6 5E6]);
- c.Label.String = 'Stress (Pa)';
- %xlabel('z (pixels)'); ylabel('y (pixels)'); axis square;
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement