Skull Shape Analysis

As a coding enthusiast, BE100 – Problem-Solving in Biological Sciences and Engineering was the class I enjoyed most during this fall quarter. The course focused on applying MATLAB to real biomedical problems, such as writing code to convert DNA sequences into amino acid chains. Through these assignments, I learned how to plot data, perform operations with matrices, and analyze biological information. All these skills eventually came together in our final project, where we basically…compare skulls. It’s obviously more than that … read on! This was technically a group project as we were supposed to write the code on our own but we could get assisted from the people at our table if we happened to be stuck. For that reason, I will acknowledge that help and use the pronouns “we/our/us” to describe the process.

What is the code doing exactly?

Overall, our code collects the skull contour, nose, and eye points from each of the given skull images. Following this pattern:

Using the points we collected, we plot the original, translated, scaled, and rotated coordinates for the skulls to then determine the smallest/largest Procrustes distance between each of the skulls and Skull One to reveal which skull is most/least like Skull One.

1: This was accomplished by using various individual functions. The Skulls_GetSkullPointsFromImage.m file enables the user to select a skull image and click on it 18 times to outline the skull contour, nose, and both eyes.

clear
clc

% Open a dialog by which the user can select the name of a jpg file that is intended to be an image of a skull.
% Store the full pathname for the file in a variable called jpg_file.
[fileName, filePath] = uigetfile('*.jpg', 'Select a file to read ...');
jpg_file = fullfile(filePath, fileName);


% Close all existing figures and start a new figure.
close all

% In the left half of a 1 x 2 subplot, show the image in the file: Skulls_18_points.jpg

figure
subplot(1, 2, 1)
hold on % Not needed for these subplots, but needed in general
data = imread("Skulls_18_points.jpg"); 
imshow(data)
title('Please Select These 18 Points ->')

% Read the image in the file whose name the user previoulsy selected.
% Store the result in a variable named something like  skull_matrix.
skull_matrix = imread(jpg_file);

% In the right half of a 1 x 2 subplot, show your skull_matrix.
subplot(1, 2, 2)
hold on % IS needed here, since two plots on the same axes
imshow(skull_matrix)
title('Your Selected Skull')

% In a loop that runs 18 times:
%   -- Allow the user to select a SINGLE point from the right image.
%   -- Display a red filled circle at the selected point (as feedback for the user).
% Hopefully the user will select 18 points in the right image as suggested by the Skulls_18_points.jpg image in the left image. As the loop runs, store the selected points in a matrix with 2 columns (for the x-coordinates and y-coordinates, respectively). So, after the user finishes selecting 18 points, you should have a variable that is an 18 x 2 matrix specifying the selected points in column vectors.
userPoints = [];
for k=1:18
    [userX, userY] = ginput(1);
    plot(userX, userY, 'or')
    userPoints = [userPoints; userX, userY];
end 

% Create a new variable that is the same as the filename that the user selected previously in this script, but ending in .xlxs instead of .jpg.
excel_filename = replace(jpg_file, 'jpg', 'xlsx');
dataRange = "A1:B18";

% Write the matrix of 18 points into an Excel file specified by your excel_filename variable.
writematrix(userPoints, excel_filename, 'Sheet',1,'Range', dataRange)

2: The Skulls_PlotPolygon.m function plots the given matrix (coordinates for a selected skull) in the given color and line width.

function Skulls_PlotPolygon(polygon_matrix, color, line_width)

% What comes in:
%  - An N x 2 matrix where each row specifies the X and Y of a point
%       (and hence the matrix specifies the coordinates of a polygon).
%  - The color for the lines in the plot of the given polygon.
%  - The line width for the lines in the plot of the given polygon.
% What goes out:
%  - Nothing (i.e. there is no return variable)
% Side effects:
%  - Draws (that is, plots) the specified polygon as a CLOSED shape.

initial_hold_state = ishold;
hold on 

% Create a variable for the SIDES (lines) of the CLOSED polygon
% by appending the first point of the polygon to the end of the
% polygon_matrix.  That way, the subsequent  plot  statement will
% include a line from the last point of the polygon back to the
% first point of the polygon, thus drawing it as a CLOSED shape.
sides = [polygon_matrix; polygon_matrix(1, :)];


% Create variables  xs  ys  for the x-coordinates (column 1)
% and y-coordintes (column 2) of the polygon_matrix.
xs = sides(:,1);
ys = sides(:,2);

% Plot the xs and ys, using the specified Color and LineWidth
plot(xs, ys, "Color", color, "LineWidth",line_width)

if initial_hold_state
    hold on
else
    hold off
end

end

3: With the existing plotting function, we use Skulls_PlotPolygon.m multiple times in Skulls_PlotSkull.m to create a brief sketch of the skulls.

function Skulls_PlotSkull(skull_matrix, color, line_width)

% What comes in:
%  - An 18 x 2 matrix where each row specifies the X and Y of a point
%      and the rows of the matrix specify the anatomical markers
%      for 4 features of a skull, as follows:
%        - Rows 1 to 7 define the skull contour
%        - Rows 8 to 10 define the nose
%        - Rows 11 to 14 define the right eye
%        - Rows 15 to 18 define the left eye
%  - The color for the lines in the plot of the given skull features.
%  - The line width for the lines in the plot of the given skull features.
% What goes out:
%  - Nothing (i.e. there is no return variable)
% Side effects:
%  - Draws (that is, plots) each of the 4 skull features
%     as a CLOSED shape, with the given color and line_width.

% HINT: call  Skulls_PlotPolygon.m   to plot each of the 4 skull features.

% TODO:
% Extract four N x 2 matrices from the given 18 x 2 skull matrix,
% defining four variables for them:
%   skull_contour = rows 1 to 7 of skull_matrix
%   nose =          rows 8 to 10 of skull_matrix
%   right_eye =     rows 11 to 14 of the skull matrix
%   left_eye =      rows 15 to 18 of the skull matrix
skull_contour = skull_matrix(1:7, :);
nose = skull_matrix(8:10, :);
right_eye = skull_matrix(11:14, :);
left_eye = skull_matrix(15:18, :);

% TODO: Plot the skull_contour as a closed polygon, with the
% given color and line_width.  (Call Skulls_PlotPolygon to do this.)
Skulls_PlotPolygon(skull_contour, color, line_width);

% TODO: Plot the nose as a closed polygon, with the
% given color and line_width.  (Call Skulls_PlotPolygon to do this.)
Skulls_PlotPolygon(nose, color, line_width);

% TODO: Plot the right_eye as a closed polygon, with the
% given color and line_width.  (Call Skulls_PlotPolygon to do this.)
Skulls_PlotPolygon(right_eye, color, line_width);

% TODO: Plot the left_eye as a closed polygon, with the
% given color and line_width.  (Call Skulls_PlotPolygon to do this.)
Skulls_PlotPolygon(left_eye, color, line_width);

% Set the axis direction and scale appropriately for converting
% from Excel's coordinate system to Matlab's coordinate system.
% Simply use the following two statements (just un-comment them):
axis ij;
axis equal;

4: After that, Skulls_TranslatePolygon.m translates the selected skull coordinates so that both images are at the origin.

function translated_polygon = Skulls_TranslatePolygon(original_polygon)

% Given a 2 x N matrix that represents N points
% (and hence a polygon), return a 2 x N matrix where each point
% of the given original_polygon is translated by an amount that
% makes the centroid of the returned polygon be at (0, 0)

centroid = mean(original_polygon);
translated_polygon = original_polygon - centroid;

end

5: And Skulls_ScalePolygon.m scales the selected skull coordinates so that both images have the same size for easier comparison.

function scaled_polygon = Skulls_ScalePolygon(original_polygon)

% Given a 2 x N matrix that represents N points
% (and hence a polygon), return a 2 x N matrix where each point
% of the given original_polygon is scaled by an amount that
% makes the size of the returned polygon be 1.
% The original_polygon is assumed to be already translated
% so that its centroid is at (0, 0)

dists = sqrt(original_polygon(:,1).^2 + original_polygon(:,2).^2);
size = sum(dists);

scaled_polygon = original_polygon * (1/size);

end

6: Then, Skulls_BestRotatedPolygon.m calculates the smallest Procrustes distance between two given polygons by rotating one degree at a time and calculates a series of Procrustes distances (Skull One and whatever other skull the user selects in Skulls_CompareSkulls).

function best_rotated_polygon = Skulls_BestRotatedPolygon(reference_polygon, polygon_to_rotate)

% Inside a loop that goes 360 times:
%   1. Rotate the  polygon_to_rotate  1 degree (use MyRotateFunction)
%   2. Compute the Procrustes Distance between the  reference_polygon
%        and the just-rotated polygon_to_rotate (use ProcrustesDistance)
%   3. Include code that finds the rotated polygon that has
%        the smallest Proscrustes Distance.

smallest_distance = inf;
for k=1:360
    rotated_polygon = MyRotateFunction(polygon_to_rotate, k); 
    procrustes_distance = ProcrustesDistance(reference_polygon, rotated_polygon);
    if procrustes_distance < smallest_distance
        smallest_distance = procrustes_distance;
        best_rotated_polygon = rotated_polygon; 
    end 
end 

end

With all the needed data generated, we collected the Procrustes distance between each of the skulls and Skull One to determine which skulls were the most similar/different.

While making this project, we learned and applied our knowledge of reading/writing Excel files, reading/showing images, and plotting. We had to read and show each of the given skull image files and then write an Excel file to store coordinates. To plot these coordinates (and their translated, scaled, and rotated counterparts), we had to read the Excel files specified by the user. More importantly, it shows that eyeballing or judging with intuition can be wrong; and to be as accurate as possible, scientific and complicated data analysis should assist in generating the right outcome. These are all skills that we learned separately through reading sample code and doing individual assignments, but this project helped us see how these three important skills can be integrated into each other. Developing these skills of making various functions that will then be integrated into a larger function will be helpful in our future BE classes, where we might be asked to analyze lab data.

An extra feature we added to the code was expanding the subplot to include 6 sections instead of 4. The two additional sections at the beginning show the two skulls being compared in the subsequent plots. We decided to add this feature because when we ran the code, we kept looking back at the image for the corresponding skull to see if the provided Procrustes distance looked reasonable. To make the program more user friendly, we added these images so the user would not have to flip back to their files to look at the corresponding skulls when running the code. Everything required for this addition is included in the image of code below from the Skulls_CompareSkulls file and produces the below image.

Skulls_CompareSkulls.m

% Before running this script, you should run your
%   Skulls_GetSkullPointsFromImage.m
% script at least twice,
% thus creating Excel files for at least two of the skull images.

clear
clc

% Open a dialog by which the user can select the name of an Excel file
% that is intended to contain the 18 points from an image of a skull.
% Store the full pathname for the file in a variable called something like

[fileName, filePath] = uigetfile('*.xlsx', 'Select a file to read ...');
skull_1_file = fullfile(filePath, fileName);

% Repeat the previous step, that is, again open a dialog
% by which the user can select the name of an Excel file
% that is intended to contain the 18 points from an image of a skull.
% Store the full pathname for the file in a variable called something like
%   skull_2_file.
[fileName, filePath] = uigetfile('*.xlsx', 'Select a file to read ...');
skull_2_file = fullfile(filePath, fileName);

% Read the two Excel files into variables something like
original_skull_1_matrix = readmatrix(skull_1_file);
original_skull_2_matrix = readmatrix(skull_2_file);

% Close all existing figures and start a new figure.
close all

% In the upper-left of a 2 x 2 subplot with title 'Original skulls',
% call  Skulls_PlotSkull  twice to plot:
%   - original_skull_1_matrix in blue
%   - original_skull_2_matrix in red
% (Choose whatever line widths you think look good.)
figure

%This code below is part of my additional piece of the project:------------ 
image_1_filename = replace(skull_1_file, 'xlsx','jpg');
image_2_filename = replace(skull_2_file, 'xlsx', 'jpg'); 

subplot(2, 3, 1)
hold on 
dataOne = imread(image_1_filename); %Where do we get this file!!!?
imshow(dataOne)
title('Skull We are Comparing (blue)')

subplot(2, 3, 2)
hold on 
dataTwo = imread(image_2_filename); %Where do we get this file!!!?
imshow(dataTwo)
title('Choosen Skull (red)')
%--------------------------------------------------------------------------

subplot(2, 3, 3)
hold on % Not needed for these subplots, but needed in general
Skulls_PlotSkull(original_skull_1_matrix, 'b', 3)
Skulls_PlotSkull(original_skull_2_matrix, 'r', 3)
title('Original Skulls')

% Call  Skulls_TranslatePolygon   twice to get two new matrices,
% each of which translates one of the original matrices so that
% its centroid is moved to (0, 0),
% storing them in variables named something like:
translated_skull_1_matrix = Skulls_TranslatePolygon(original_skull_1_matrix);
translated_skull_2_matrix = Skulls_TranslatePolygon(original_skull_2_matrix);

% In the upper-right of a 2 x 2 subplot with title 'Translated skulls',
% call  Skulls_PlotSkull  twice to plot:
%   - translated_skull_1_matrix in blue
%   - translated_skull_2_matrix in red
% (Choose whatever line widths you think look good.)
subplot(2, 3, 4)
hold on 
Skulls_PlotSkull(translated_skull_1_matrix, 'b', 3)
Skulls_PlotSkull(translated_skull_2_matrix, 'r', 3)
title('Translated Skulls')

% Call  Skulls_ScalePolygon   twice to get two new matrices,
% each of which scales one of the translated matrices so that
% its size is 1,
% storing them in variables named something like:
scaled_skull_1_matrix = Skulls_ScalePolygon(translated_skull_1_matrix);
scaled_skull_2_matrix = Skulls_ScalePolygon(translated_skull_2_matrix);

% In the lower-left of a 2 x 2 subplot with title 'Scaled skulls',
% call  Skulls_PlotSkull  twice to plot:
%   - scaled_skull_1_matrix in blue
%   - scaled_skull_2_matrix in red
% (Choose whatever line widths you think look good.)
subplot(2, 3, 5)
hold on 
Skulls_PlotSkull(scaled_skull_1_matrix, 'b', 3)
Skulls_PlotSkull(scaled_skull_2_matrix, 'r', 3)
title('Scaled Skulls')

% Call  Skulls_BestRotatedPolygon,
%   sending it the  scaled_skull_matrix_1  as the reference skull
%   and the  scaled_skull_matrix_2  as the skull to rotate
% storing the result in a variable named something like:
rotated_skull_2_matrix = Skulls_BestRotatedPolygon(scaled_skull_1_matrix, scaled_skull_2_matrix);


% In the lower-right of a 2 x 2 subplot with title 'Rotated skulls',
% call  Skulls_PlotSkull  twice to plot:
%   - scaled_skull_1_matrix in blue
%   - rotated_skull_2_matrix in red
% (Choose whatever line widths you think look good.)
subplot(2, 3, 6)
hold on 
Skulls_PlotSkull(scaled_skull_1_matrix, 'b', 3)
Skulls_PlotSkull(rotated_skull_2_matrix, 'r', 3)
title('Rotated Skulls')