Analysis of Scoliosis X-Rays and Spinal Deviation

By: Euressa Cosmiano, Trevor Pearce, Ambrocio Sanchez, Hung Hoang

This MATLAB code was built to analyze patient x-rays, specifically those who suffer from scoliosis.

Scoliosis is a condition where the spine contains a sideways curve; the severity of this disease can range from mild deviation to extreme and disabling pain. In order to determine the severity of each patient's case, we carefully examine the curve and determine which course of treatment is best for the patient. Since each case is different, we measure the angle at which the spine is curving; this angle is also known as Cobb's angle. From there, we can address the patient with their treatment options.

In this code, we upload the specific x-ray, create a graphical representation of their spinal deviation, take a measurement of Cobb's angle, and create "if-else" statements that directly tell the user what they can do to improve this condition.

Step 1: Selecting and Uploading X-ray Images

Before starting the actual code, we began looking for clear x-ray images of patients with scoliosis. If we could not find the specific images we were looking for, then we would have had to start this assignment with a new idea. The images we were searching for had to meet specific criteria to be suitable for this project. We began looking in Google images under specific title "Scoliosis X-ray Images". We had a variety of images to choose from, so the image criteria became more strict. Each x-ray had to be a "clean copy"; there should be no evidence of previous modifications, marks, or measurements from previous doctors or medical staff. The x-rays were also required to already be black and white, similarly to a traditional x-ray. Although present day technology allows for more involved and informative x-rays, our group wanted to stay consistent throughout each image. From our narrowed selection, we chose x-rays that extended from the top of the shoulders to the hips. Next, we ensured that each vertebrae could be seen from the x-ray. We then gathered these images into a folder and batch uploaded this into MATLAB. The code begins calling upon the images within the folder and allows the user to select which image they would like to see. The image then appears in a figure along with instructions to determine Cobb's angle. ___________________________________________________________________________

%% Get Image and Display Image

GetDir = uigetfile('*.jpg');

SourceImg = imread(GetDir);

figure(1);

subplot(1, 2, 1);

imshow(SourceImg);

title('Scoliosis X-Ray');

waitfor(msgbox({'1. Expand figure window to enlarge image'; '2. Click on center of each vertebrae'; '3. Hit enter to display plot'}));

hold on

Step 2: Identifying Vertebrae and Showing Deviations

Once the image is uploaded, the user selects the center of each vertebrae by using the cursor by using the function 'ginput'. After identifying each vertebrae, the user will press 'Enter' on the keyboard. The x-ray will now contain red circles where the user clicked and a graphical representation of the curved spine will appear in the same figure, adjacent to the x-ray. In the graphical representation, there will be a "zero line" that serves as a normal spine alignment.

The user can now see where the spine deviates from the zero line, how large the curve is, and where the spine returns to normal alignment.

________________________________________________________________________________________

%% Identify Vertebrae and Show Deviations

[x, y] = ginput;

plot(x, y, 'ro', 'LineWidth', 2, 'MarkerSize', 15);

points = [x, y];

points = round(points);

xCenter = points(find(y==max(y)), 1);

xNew = x - xCenter;

subplot(1, 2, 2);

stem(y, xNew);

camroll(-90);

yticks([]);

xticks([]);

xlim([min(y)-50 max(y)+50]);

if max(abs(xNew)) >= min(abs(xNew))

ylim([-max(abs(xNew))-50 max(abs(xNew))+50]);

else

ylim([-min(abs(xNew))-50 min(abs(xNew))+50]);

end

title('Deviations From Normal Spine');

pause(2.0);

Step 3: Display Polynomial Fit With Reference Line

This section of the codes computes automatically. From the graphical representation, the code will analyze the data and take the best fit line, using the function 'polyfit'. This new line will be displayed in a figure along with the reference line for normal spine alignment.

________________________________________________________________________________________

%% Display Polynomial Fit w/ Reference Line

% Polynomial

imgSize = size(SourceImg);

p = polyfit(y, xNew, 3);

x = linspace(-50, imgSize(1), imgSize(1)+50);

f1 = polyval(p, x);

figure(2);

plot(f1, 'r', 'LineWidth', 10);

yticks([]);

xticks([]);

title('Severity of Spinal Curvature');

xlim([min(y)-50 max(y)+50]);

if max(abs(xNew)) >= min(abs(xNew))

ylim([-max(abs(xNew))-500 max(abs(xNew))+500]);

else

ylim([-min(abs(xNew))-500 min(abs(xNew))+500]);

end

hold on

% Zero line

f2 = polyval(0, x);

plot(f2, 'r');

hold on

Step 4: Find Critical Points and Normal Lines

On the same figure of the reference line, the code detects the critical points, the points where the 'polyfit' line crosses the zero line, and calculates the slopes to these lines. The code is not designed to plot the slope lines, therefore they do not show up in the figure. The next step is to take the negative reciprocal of the slope to create normal lines; these lines are then displayed onto the figure.

The intersection of the normal lines is accounted for and the plotted lines stop at their intersection, which is marked by a point.

________________________________________________________________________________________

%% Find Critical Points and Normal Lines

% Critical points

CP = sort(roots(p), 'ascend');

plot(CP(1)+50, 0, 'bo', 'LineWidth', 2, 'MarkerSize', 10); % top critical point

hold on

plot(CP(2)+50, 0, 'mo', 'LineWidth', 2, 'MarkerSize', 10); % middle critical point

hold on

plot(CP(3)+50, 0, 'ro', 'LineWidth', 2, 'MarkerSize', 10); % bottom critical point

% Slope of polynomial at critical points

x2 = linspace(CP(1)+45, CP(1)+55, 2);

f3 = polyval(p, x2);

slope1 = 10/(f3(2)-f3(1));

x3 = linspace(CP(2)+45, CP(2)+55, 2);

f4 = polyval(p, x3);

slope2 = 10/(f4(2)-f4(1));

x4 = linspace(CP(3)+45, CP(3)+55, 2);

f5 = polyval(p, x4);

slope3 = 10/(f5(2)-f5(1));

% Top intersection point and normals

m1 = -(slope1)^-1;

m2 = -(slope2)^-1;

m3 = -(slope3)^-1;

b1 = CP(1)+50;

b2 = CP(2)+50;

b3 = CP(3)+50;

A = [m2, 1; m1, 1];

b = [b2; b1];

topIP = A\b; % intersection point of normals

plot(topIP(2), -topIP(1), 'bo', 'LineWidth', 2, 'MarkerSize', 5);

hold on

if topIP(1) < 0

x5 = linspace(0, -topIP(1), -topIP(1));

else

x5 = linspace(-topIP(1), 0, topIP(1));

end

norm1 = m1*x5 + b1;

plot(norm1, x5, 'b'); % bottom line

hold on

norm2 = m2*x5 + b2;

plot(norm2, x5, 'b'); % top line

hold on

% Bottom intersection point and normals

A = [m3, 1; m2, 1];

b = [b3; b2];

botIP = A\b; % intersection point of normals

plot(botIP(2), -botIP(1), 'ro', 'LineWidth', 2, 'MarkerSize', 5);

hold on

if botIP(1) < 0

x6 = linspace(0, -botIP(1), -botIP(1));

else

x6 = linspace(-botIP(1), 0, botIP(1));

end

norm2 = m2*x6 + b2;

plot(norm2, x6, 'r'); % bottom line

hold on

norm3 = m3*x6 + b3;

plot(norm3, x6, 'r'); % top line

hold on

camroll(-90);

Step 5: Find Side Lengths of Triangle and Calculate Angle of Severity

This is the final step before calculating Cobb's Angle/ Angle of Severity.

The code analyzes the data of the critical points and the normal lines. It first takes the measurement of each normal line until its intersection with the other line. Then by using the Law of Cosines, we can calculate the two angles closest to the zero line. From here, the code will calculate the missing angle, which is the Angle of Severity. Although the angles closest to the zero line, the "Lower Angle" and "Upper Angle", are not displayed in the figure, Cobb's angle will be shown next to the triangle.

________________________________________________________________________________________

%% Find Side Lengths of Triangle and Calculate Angle of Severity

% Upper Angle

side1 = CP(2) - CP(1);

side2 = sqrt((b1 - topIP(2))^2 + (-topIP(1))^2);

side3 = sqrt((b2 - topIP(2))^2 + (-topIP(1))^2);

radians = acos(((side2)^2 + (side3)^2 - (side1)^2)/(2 * side2 * side3)); % Law of Cosines

degreesTop = radians * (180/pi);

degreesTop = round(degreesTop);

fprintf('The severity of upper spinal curvature for this patient is approximately %.f degrees', degreesTop)

str = ([num2str(degreesTop), '\circ']);

if topIP(1) < 0

text(topIP(2), -topIP(1)-60, str, 'Interpreter', 'tex', 'FontSize', 18);

else

text(topIP(2), -topIP(1)+60, str, 'Interpreter', 'tex', 'FontSize', 18);

end

% Lower Angle

side1 = CP(3) - CP(2);

side2 = sqrt((b2 - botIP(2))^2 + (-botIP(1))^2);

side3 = sqrt((b3 - botIP(2))^2 + (-botIP(1))^2);

radians = acos(((side2)^2 + (side3)^2 - (side1)^2)/(2 * side2 * side3)); % Law of Cosines

degreesBottom = radians * (180/pi);

degreesBottom = round(degreesBottom);

fprintf('\n The severity of lower spinal curvature for this patient is approximately %.f degrees', degreesBottom)

str = ([num2str(degreesBottom), '\circ']);

if botIP(1) < 0

text(botIP(2), -botIP(1)-60, str, 'Interpreter', 'tex', 'FontSize', 18);

else

text(botIP(2), -botIP(1)+60, str, 'Interpreter', 'tex', 'FontSize', 18);

end

Step 6: Diagnosis and Treatment Options Based on the More Severe Angle

This section of the code contains "if-else" statements. Through research of the condition and possible treatment plans, we organized this information according to the severity of Cobb's Angle.

Although there are general ways of dealing with this condition, we included the most realistic and advantageous treatments for varying ranges of spinal deviation. After the calculation of Cobb's Angle, the code will recognize the range in which this patient's spinal deformation belongs and return treatment plans to the patient.

________________________________________________________________________________________

%% Diagnosis and Treatment Options Based on the More Severe Angle

if degreesTop >= degreesBottom

degrees = degreesTop;

else

degrees = degreesBottom;

end

if degrees > 40

fprintf('Please schedule the patient for a surgery as soon as possible.\nThere are multiple types of surgery available:\n1. Spinal fusion and metal rod insertion.\n2. Posterior approach.\n3. Anterior approach.\n4. Decompressive laminectomy.\n5. Minimally invasive surgery.\n')

fprintf('Please keep in mind that if the patient is not an adult,\nthey may require revision surgery in 20 - 30 years. In the case\nthat newer treatment is avaialable.\n')

elseif degrees > 20 && degrees <= 40

fprintf('Please recommend that the patient wear a brace for 16 - 23 hours a day.\nIf the patient is a child or teen, the recommended duration of the brace is until\nage of maturity.\n If the patient is an adult, wearing a brace might not help\nmuch, so surgery might be the better option.\n')

elseif degrees > 10 && degrees <= 20

fprintf('Please recommend that the patient schedules a regular check-in every 2 years.\nProgression of scoliosis may worsen over time, but at the moment the patient\ndoes not require any treatment.\n')

end