In my previous post I detailed methods for conducting a binary image segmentation of the corpus callosum. I used three different image modalities (T1-weighted, T2-weighted and generalized fractional anisotropy) and a simple k-Nearest Neighbors (kNN) model from sklearn. We also covered some basic normalization that is important for kNN models.
Today’s goal is to demonstrate how some knowledge of your dataset can help you derive informative features.
In our dataset, we need to deal with some misclassified tissue from my last post on the corpus callosum. As a reminder of our results, lets load our data from last time and visualize the results.
from nilearn import plotting
import matplotlib.pyplot as plt
import pandas as pd
from sklearn import preprocessing
from sklearn.utils import resample
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split
import numpy as np
import glob
import warnings
warnings.filterwarnings('ignore')First, a reminder of our MRI datatypes.
# create a figure with subplots
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15,4));
# plot the T1-weighted image
plotting.plot_anat("data/subject1/T1.nii.gz",
draw_cross=False,
title='T1-weighted',
axes=axes[0],
annotate=False)
# plot the T2-weighted image
plotting.plot_anat("data/subject1/T2.nii.gz",
draw_cross=False,
title='T2-weighted',
axes=axes[1],
annotate=False)
# plot the GFA image
plotting.plot_anat("data/subject1/GFA.nii.gz",
draw_cross=False,
title='GFA',
axes=axes[2],
annotate=False)<nilearn.plotting.displays.OrthoSlicer at 0x11a691630>
Again, I’ve cheated a bit here and I’ve previously done some extensive preprocessing to this dataset. Stay tuned for a post detailing that process.
As before, let’s load the files, normalize that data, deal with our class imbalance and train a series of models to predict left out data. For slightly more information please revisit the previous post.
# load all data files
df = pd.DataFrame()
all_files = glob.glob('data/*.limitedData.csv')
for ff in all_files:
data = pd.read_csv(ff)
df = pd.concat([df, data], ignore_index=True)
# separate response and predictor data
Y = pd.Series(df["response"])
X = df.drop(labels='response', axis=1)
# scale
colnames = list(X)
scaler = preprocessing.RobustScaler()
scaled_X = scaler.fit_transform(X)
scaled_X = pd.DataFrame(scaled_X, columns=colnames)# create new response and predictor variables
X_class0 = scaled_X[Y == 0]
X_class1 = scaled_X[Y == 1]
# downsample predictors to equal the minority class (1's)
X_downSampled = resample(X_class0,
replace = True,
n_samples = len(scaled_X[Y == 1]),
random_state = 99)
# remerge the data
Y_matched = pd.concat([pd.DataFrame(0, index = np.arange(len(X_downSampled),1), columns=['Response']), Y[Y == 1]], ignore_index=True, axis = 0)
X_matched = pd.concat([X_downSampled, scaled_X[Y == 1]], ignore_index=True)# create a subject identifier to allow for leave-one-subject-out loop
subject_identifier = np.zeros(len(scaled_X));
start = 0
stop = 4500
for ss in range(1,len(all_files)+1):
subject_identifier[start:stop] = (np.ones(4500, dtype = int) * ss)
start = start + 4500
stop = stop + 4500# loop over subjects, leaving one out each time
for ss in range(1,len(all_files)+1):
# split, train, test
X_train, X_test, Y_train, Y_test = train_test_split(scaled_X[subject_identifier!=ss],Y[subject_identifier!=ss], test_size = 0.2, random_state = 1)
knn = KNeighborsClassifier(n_neighbors = 5)
mdl = knn.fit(X_train, Y_train)
# for 3 example subjects ...
if ss < 4:
preds = mdl.predict(scaled_X[subject_identifier==ss])
preds = np.reshape(preds,(50,90))
test = Y[subject_identifier==ss].values.reshape(50,90)
# reorient the images
reorientedPreds = np.fliplr(preds)
reorientedPreds = np.flipud(reorientedPreds)
reorientedTest = np.fliplr(test)
reorientedTest = np.flipud(reorientedTest)
# plot the actual and predicted
fig = plt.figure(figsize=(10,3))
a = fig.add_subplot(1, 2, 1)
imgplot = plt.imshow(reorientedPreds)
a.set_title('Predicted')
a = fig.add_subplot(1, 2, 2)
imgplot = plt.imshow(reorientedTest)
a.set_title('Actual')
plt.show()
As we can see from our predicted class membership, we are doing well with the body of the corpus callosum but are misclassifying voxels near the edges and in a nearby structure called the fornix. This occurs, in part, because of the limitations of our current predictors. T1-weighted and T2-weighted images are roughly anticorrelated, and do a good job of distinguishing gray matter (primarily the neurons) from the white matter (primarily the axons). Note that this is a significant oversimplification of the number and variety of cells in the human brain, but serves to illustrate my point.
GFA is derived from diffusion imaging, which samples the microscopic movement of water molecules in the brain. Larger GFA values indicate that the movement has a principle direction (e.g. along the direction of a bundle of axons) while lower GFA values indicate more random movement (e.g. within cerebrospinal fluid). The information from the GFA data also aids in the segmentation of gray matter, white matter and CSF as they each have different water diffusion properties.
Now one thing that GFA doesn’t do is distinguish between white matter containing axons that have different principle directions. For example, axons travelling rostrally (towards the front of your head) will have similiar GFA values to voxels containing axons travelling dorsally (towards the top of your head), despite their different principal directions. So when we use GFA in our machine learning models there is no way for the model to learn to distinguish between different axon bundles that might be present in the image together!
Luckily we can use our knowledge of anatomy to develop some features that might help distinguish the corpus callosum from the fornix!
Deriving additional informative features
From our super advanced knowledge of anatomy, we know that the axons of the corpus callosum primarily travel left-to-right across the brain. This is especially true at the middle of the brain where our data is sampled from. We also know that the fibers of the fornix primarily travel first medially (towards the middle of your head) and then laterally (towards the sides of your head) as the proceed ventral-to-dorsal - but their principle direction is roughly orthogonal to the corpus callosum axons. If we could come up with a feature, or set of features, that captures this information we might be able to harness that information for our model.
This is where a technique called diffusion tensor imaging (DTI) can come in handy. For simplicity’s sake I won’t go into the details of DTI - in fact the DTI analysis is how I initially calculated the GFA values. What DTI allows us to do, however, is estimate the principle direction of the axons within a voxel according to two vectors phi and theta in spherical polar coordinates (see this overview for more information). We should be able to use these two values to come up with the features we need.
Let’s look at the data!
# create a figure with subplots
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15,4));
# plot the GFA image
plotting.plot_anat("data/subject1/GFA.nii.gz",
draw_cross=False,
title='GFA',
axes=axes[0],
annotate=False)
# plot the phi image
plotting.plot_anat("data/subject1/phi.nii.gz",
draw_cross=False,
title='phi',
axes=axes[1],
annotate=False)
# plot the theta image
plotting.plot_anat("data/subject1/theta.nii.gz",
draw_cross=False,
title='theta',
axes=axes[2],
annotate=False)
## copy the phi and theta images for subject 1<nilearn.plotting.displays.OrthoSlicer at 0x11a2bab70>
The phi and theta images look like a lot of noise, but if you look closely you can see that the corpus callosum is visible as a brighter structure in the phi data and a darker structure in the theta data. The converse is true for the fornix.
As I said before, the phi and theta images represent the principle fiber direction for axons contained in a given voxel. Originally these values are in radians, so to derive a ‘nice’ metric from them I first took the cosine of the value for each voxel to scale it between 0 and 1, and then, because the tool I used to calculate the cosine of the image values maintained the sign, I also took the absolute value. This had the benefit of removing negative values which were really just the opposite angle in 180 degree space (e.g. 90 degress and 270 degrees are really the same thing for our purposes).
So we’ve successfully derived new features from our data that behave nicely and would seem to give some useful information for the model to take advantage of.
Classification
Let’s try using these as features but simply adding them to our feature set. Again, I’ve previously collated these data points for each subject in our dataset.
# load all data files
df = pd.DataFrame()
all_files = glob.glob('data/*.allData.csv')
for ff in all_files:
data = pd.read_csv(ff)
df = pd.concat([df, data], ignore_index=True)
# separate response and predictor data
Y = pd.Series(df["response"])
X = df.drop(labels='response', axis=1)
# scale
colnames = list(X)
scaler = preprocessing.RobustScaler()
scaled_X = scaler.fit_transform(X)
scaled_X = pd.DataFrame(scaled_X, columns=colnames)# create new response and predictor variables
X_class0 = scaled_X[Y == 0]
X_class1 = scaled_X[Y == 1]
# downsample predictors to equal the minority class (1's)
X_downSampled = resample(X_class0,
replace = True,
n_samples = len(scaled_X[Y == 1]),
random_state = 99)
# remerge the data
Y_matched = pd.concat([pd.DataFrame(0, index = np.arange(len(X_downSampled),1), columns=['Response']), Y[Y == 1]], ignore_index=True, axis = 0)
X_matched = pd.concat([X_downSampled, scaled_X[Y == 1]], ignore_index=True)# loop over subjects, leaving one out each time
for ss in range(1,len(all_files)+1):
# split, train, test
X_train, X_test, Y_train, Y_test = train_test_split(scaled_X[subject_identifier!=ss],Y[subject_identifier!=ss], test_size = 0.2, random_state = 1)
knn = KNeighborsClassifier(n_neighbors = 5)
mdl = knn.fit(X_train, Y_train)
# plot 3 example left out subjects
if ss < 4:
predsAllData = mdl.predict(scaled_X[subject_identifier==ss])
predsAllData = np.reshape(predsAllData,(50,90))
test = Y[subject_identifier==ss].values.reshape(50,90)
# reorient the images
reorientedPredsAllData = np.fliplr(predsAllData)
reorientedPredsAllData = np.flipud(reorientedPredsAllData)
reorientedTest = np.fliplr(test)
reorientedTest = np.flipud(reorientedTest)
# plot the actual and predicted
fig = plt.figure(figsize=(10,3))
a = fig.add_subplot(1, 2, 1)
imgplot = plt.imshow(reorientedPredsAllData)
a.set_title('Predicted, All Data')
a = fig.add_subplot(1, 2, 2)
imgplot = plt.imshow(reorientedTest)
a.set_title('Actual')
plt.show()
No more fornix! Huzzah! The information the model has incorporated from the phi and theta images were enough to remove the misclassified tissue.
So let’s recap. So far we’ve 1) loaded data and used simple, classic machine learning methods to obtain a reasonably good approximation of the corpus callosum, 2) refined our methods to include normalization and leave-one-out cross validation and 3) derived features from the dataset of aid in our classification/segmentation efforts.
There is still some work to be done here. Especially around the edges of the corpus callosum where the values are more varied. And there is also some schmutz at the bottom of the image. In my next post I will detail some basic computer vision techniques for removing these misclassified voxels before we wrap this up into a tidy function. Then maybe I’ll get around to writing that post on MRI preprocessing…