In [35]:
import os
import numpy as np # linear algebra
import scipy
from skimage.io import imread
import matplotlib.pyplot as plt
from skimage import segmentation
from skimage.filters import gaussian, median
from skimage.measure import *
from sklearn.cluster import *
from skimage.segmentation import *
from skimage.exposure import *
import matplotlib.patches as patches
from skimage.transform import *
from matplotlib.path import Path
import peakutils
from peakutils.plot import plot as pplot
import string
import csv
In [36]:
# cut machine information from image
def crop_border(image):
x1,x2 = 120,571
y1,y2 = 60,531
image = image[y1:y2,x1:x2]
return resize(image,[451,441])
In [37]:
# cut training image to match
def crop_border_ref(image):
x1,x2 = 203,644
y1,y2 = 87,538
return image[y1:y2,x1:x2]
In [38]:
IMAGE_BICEP_TRAIN = './tendon_images/Bicep_Train'
IMAGE_SUP_TRAIN = './tendon_images/Sup_Train'
IMAGE_BICEP = './tendon_images/Biceps'
IMAGE_SUP = './tendon_images/Supraspinatus'
def load_images_names_from_folder(folder):
images = []
name = []
train = []
train_name = []
for filename in os.listdir(folder):
if "ID" in filename and 'REF' not in filename:
img = imread(os.path.join(folder,filename))
images.append(crop_border(img))
name.append(filename)
# if there are training images, load them
ref = filename.replace('.jpg','')+'_REF_SEL.jpg'
if ref in os.listdir(folder):
img = imread(os.path.join(folder,ref))
train.append(crop_border_ref(img))
train_name.append(ref)
return [images,train,name,train_name]
[images,train,names,train_name] = load_images_names_from_folder(IMAGE_SUP_TRAIN)
In [39]:
# plot test image
%matplotlib inline
fig, (ax1,ax2) = plt.subplots(1,2, figsize = (5,5))
ax1.imshow(images[5])
ax1.set_title('Ultrasound Image')
if train:
ax2.imshow(train[2])
ax2.set_title('Reference image')
In [40]:
# extract x value to find where the ROI should start
# this looks at the very top of the image to find the imprint of two metal bars on the skin
# these are shown as peaks in the rows and we can localize an x coordinate for the region from these.
# Bicep tendon has bars on the left, sup on the right
# for the Bicep the ROI is then a set distance from the bars, for the sup it varies between patients
def set_x_ROI(image,offset):
image = gaussian(image,sigma=1)
y = []
for i in range(image.shape[1]):
y.insert(i,np.mean(image[0:10,i])) # just use top 10 rows
x = np.linspace(0, image.shape[1], image.shape[1])
y = np.array(y)
indexes = np.array(peakutils.indexes(y, thres=0.8, min_dist=20)) # there should be two peaks it returns
print indexes
x = ((indexes[-2] - indexes[-1]))
x = indexes[-2] - x - offset
#print indexes
return x
In [41]:
# extract initial ROI by sliding window down to find maximum intensity in thresholded image
# we also slide n pixels around the init extracted ROI since the ROI can change for the sup tendon
def find_initial_ROI(image,init_x):
maximum = 0
max_x = init_x
max_y = 50
(x1,y1) = (init_x,50) # if there is a brighter region at top from skin then ignore
length = 93 # from training
width = 100
(x2,y2) = (x1+length,y1+width)
#print
#print init_x
for i in range ( init_x-50, init_x+50) :
#print i
for j in range (image.shape[1] - width):
x = i
y = y1 + j
(x2,y2) = (x+length,y+width)
temp = np.sum(image[y:y2, x:x2])
if temp > maximum:
maximum = temp
max_x = x
max_y = y
#print max_x
return max_x,max_y
In [42]:
# fit contour on region based on window on maximum
def fit_contour(image,max_x,max_y,length,width):
s = np.linspace(0, 2*np.pi, 400)
x = (max_x + length/2) + length*np.cos(s)
y = (max_y + width/2) + (width/4)*np.sin(s)
init = np.array([x, y]).T
snake = active_contour(image,init, bc ='fixed', alpha=1.5, beta=1.0, w_line=1, w_edge=1, gamma=0.1)
return snake
In [43]:
# fix x coordinates of snake to be maximum length
def fix_snake(snake,x,length):
for val in snake:
if val[0] < x:
val[0] = x
if val[0] > (x+length):
val[0] = x+length
return snake
In [44]:
# based on the snake coordinates, we can extract the maximim sized
# rectangle for the region
def ROI_from_snake(snake):
s = snake
x = []
y = []
top = []
bottom = []
for i,j in s:
x.append(i)
y.append(j)
for i,j in s:
if (i >= (min(x) + 10)) and (i <= (max(x) - 10)):
if j > int(np.mean(y)):
bottom.append(int(j))
elif j < int(np.mean(y)):
top.append(int(j))
ymin = int(np.mean(y))
ymax = int(np.mean(y))
for i in range(100):
ymin = ymin - 1
if ymin in top:
break
for i in range(100):
ymax = ymax + 1
if ymax in bottom:
break
ROI = [[int(min(x)),ymin],[int(max(x)),ymax] ]
return ROI
In [45]:
# 1. Transform images by equalizing histogram and filtering
# 2. Then threshold the image to get the maximum
# 3. Extract the inital ROI
# 4. Build an intermediate thresholded image at a lower threshold
# this is to allow a better fit on the snake
# 5. Fit a snake onto the new thresholded image
# 6. Cut the snake so it fits within in the desired length
thres_val = 0.96
intermediate_thres_val = 0.92
length = 97 #1cm
width = 100
offset = 250
c = [] # converted images
snake = []
ROI_xy = []
ROI = []
rect = []
circ = []
x_vals = []
ROI_img = [] # used for F1 extraction
Precision = []
Recall = []
F1 = []
inter = [] # intermediate images to fit snake
for i in range(len(images)):
temp = images[i]
temp = equalize_hist(temp)
temp = gaussian(temp,sigma=2)
x = np.max(temp,2)
x[x<thres_val]=0
c.insert(i,x)
for i in range(len(images)):
temp = images[i]
temp = equalize_hist(temp)
temp = gaussian(temp,sigma=1)
x = np.max(temp,2)
x[x<intermediate_thres_val]=0
inter.insert(i,x)
for i in range(len(c)):
ROI_x = set_x_ROI(images[i],offset)
x_vals.insert(i,ROI_x)
[x,y] = find_initial_ROI(c[i],ROI_x)
ROI_xy.insert(i,[x,y])
contour = fit_contour(inter[i],ROI_xy[i][0],ROI_xy[i][1],length,width)
snake.insert(i,contour)
snake[i] = fix_snake(snake[i],ROI_xy[i][0],length)
region = ROI_from_snake(snake[i])
ROI.insert(i,region)
#for i in range(len(ROI)):
# r= patches.Rectangle((ROI[i][0][0],ROI[i][0][1]),length,(ROI[i][1][1] - ROI[i][0][1]),linewidth=1,edgecolor='r',facecolor='none')
# rect.insert(i,r)
for i in range(len(c)):
r= patches.Rectangle((ROI_xy[i][0],ROI_xy[i][1]),length,width,linewidth=1,edgecolor='r',facecolor='none')
rect.insert(i,r)
In [46]:
# original rectangle found
# just for plotting purpose
for i in range(len(ROI)):
r= patches.Rectangle((ROI[i][0][0],ROI[i][0][1]),length,(ROI[i][1][1] - ROI[i][0][1]),linewidth=1,edgecolor='r',facecolor='none')
rect.insert(i,r)
# this large rectanlge show the search region
large_rect = []
for i in range(len(c)):
m= plt.Rectangle((x_vals[i]-50,50),length+100,451-50,linewidth=1,edgecolor='g',fill=False)
large_rect.insert(i,m)
# plot one set of images in greater detail
%matplotlib inline
fig, ((ax1,ax2)) = plt.subplots(1,2, figsize = (14,14))
idx = 5
ax1.imshow(images[idx], cmap = 'gray')
ax1.set_title('Ultrasound Image')
ax2.imshow(c[idx], cmap = 'gist_earth')
ax2.plot(snake[idx][:, 0], snake[idx][:, 1], '-b', lw=3)
ax2.add_artist(rect[idx])
ax2.add_artist(large_rect[idx])
ax2.imshow(inter[idx], alpha=0.3, cmap='Blues')
#ax2.imshow(train[idx], alpha=0.3, cmap='Blues')
ax2.set_title('Threshold and train image with ROI and search region')
Out[46]:
In [47]:
# original rectangle found
# just for plotting purpose
for i in range(len(c)):
r= patches.Rectangle((ROI_xy[i][0],ROI_xy[i][1]),length,width,linewidth=1,edgecolor='r',facecolor='none')
rect.insert(i,r)
for i in range(len(ROI)):
r= patches.Rectangle((ROI[i][0][0],ROI[i][0][1]),length,(ROI[i][1][1] - ROI[i][0][1]),linewidth=1,edgecolor='r',facecolor='none')
rect.insert(i,r)
large_rect = []
for i in range(len(c)):
m= plt.Rectangle((x_vals[i]-50,50),length+100,451-50,linewidth=1,edgecolor='g',fill=False)
large_rect.insert(i,m)
# plot all 6 sets of images
%matplotlib inline
fig, ((ax1,ax2),(ax4,ax5),(ax7,ax8),(ax10,ax11),(ax13,ax14),(ax16,ax17)) = plt.subplots(6,2, figsize = (20,20))
# first image
ax1.imshow(images[0], cmap = 'gray')
ax2.imshow(c[0], cmap = 'gist_earth')
ax2.plot(snake[0][:, 0], snake[0][:, 1], '-b', lw=3)
ax2.add_patch(rect[0])
ax2.add_artist(large_rect[0])
# second image
ax4.imshow(images[1], cmap = 'gray')
ax5.imshow(c[1], cmap = 'gist_earth')
ax5.plot(snake[1][:, 0], snake[1][:, 1], '-b', lw=3)
ax5.add_patch(rect[1])
ax5.add_artist(large_rect[1])
# third image
ax7.imshow(images[2], cmap = 'gray')
ax8.imshow(c[2], cmap = 'gist_earth')
ax8.plot(snake[2][:, 0], snake[2][:, 1], '-b', lw=3)
ax8.add_patch(rect[2])
ax8.add_artist(large_rect[2])
# fourth image
ax10.imshow(images[3], cmap = 'gray')
ax11.imshow(c[3], cmap = 'gist_earth')
ax11.plot(snake[3][:, 0], snake[3][:, 1], '-b', lw=3)
ax11.add_patch(rect[3])
ax11.add_artist(large_rect[3])
# fifth image
ax13.imshow(images[4], cmap = 'gray')
ax14.imshow(c[4], cmap = 'gist_earth')
ax14.plot(snake[4][:, 0], snake[4][:, 1], '-b', lw=3)
ax14.add_patch(rect[4])
ax14.add_artist(large_rect[4])
# sixth image
ax16.imshow(images[5], cmap = 'gray')
ax17.imshow(c[5], cmap = 'gist_earth')
ax17.plot(snake[5][:, 0], snake[5][:, 1], '-b', lw=3)
ax17.add_patch(rect[5])
ax17.add_artist(large_rect[5])
# these show the training image. Comment out if not used
ax2.imshow(train[0],alpha=0.3)
ax5.imshow(train[1],alpha=0.3)
ax8.imshow(train[2],alpha=0.3)
ax11.imshow(train[3],alpha=0.3)
ax14.imshow(train[4],alpha=0.3)
ax17.imshow(train[5],alpha=0.3)
Out[47]:
In [48]:
# Find F1, precision, and recall scores for the images
# this is done by converting the ROI region and the training image
# to a boolean matrix and checking the overlaps.
# build image
for idx in range(len(c)):
image = c[idx]
x1 = ROI[idx][0][0]
x2 = ROI[idx][1][0]
y1 = ROI[idx][0][1]
y2 = ROI[idx][1][1]
# build rectangle path around ROI
p = Path([[x1,y1],[x2,y1],[x2,y2],[x1,y2]])
#p = Path(np.array(snake[idx]))
for i in range(image.shape[0]):
for j in range(image.shape[1]):
if (p.contains_points([[j,i]])):
image[i,j] = 255
else:
image[i,j] = 0
ROI_img.insert(idx,image)
# find f1 score
for idx in range(len(c)):
t = np.max(train[idx],2).astype(bool).ravel() # convert to bool vector
img = ROI_img[idx].astype(bool).ravel()
tp = 0.0
fn = 0.0
tn = 0.0
fp = 0.0
for i in range(len(t)):
if ((t[i] == True) and (img[i] == True)):
tp = tp +1
if ((t[i] == True) and (img[i] == False)):
fn = fn +1
if ((t[i] == False) and (img[i] == False)):
tn = tn +1
if ((t[i] == False) and (img[i] == True)):
fp = fp +1
print tp, fp, fn, tn
precision = tp/ (tp + fp)
recall = tp/(tp + fn)
Precision.insert(idx,precision)
Recall.insert(idx,recall)
f1 = 2 * (precision*recall) / (precision + recall)
F1.insert(idx,f1)
print precision, recall, f1
print
In [49]:
fig, (ax1,ax2) = plt.subplots(1,2, figsize = (10,10))
idx = 4
# first image
ax1.imshow(ROI_img[idx], cmap = 'gray')
ax1.imshow(train[idx],alpha=0.3)
ax1.set_title('Training image and extracted ROI ')
ax2.imshow(ROI_img[1], cmap = 'gray')
ax2.imshow(train[1],alpha=0.3)
ax2.set_title('Training image and extracted ROI ')
Out[49]:
In [81]:
# Save the data to a csv file
label = ['name','ROIx1','ROIy1','ROIx2','ROIy2','tendon_width','F1','Prec','Recall']
data = []
for i in range(len(names)):
x = [names[i].replace('.jpg',''),int(ROI[i][0][0]),
int(ROI[i][0][1]),int(ROI[i][1][0]),int(ROI[i][1][1]),
int(ROI[i][1][1] - ROI[i][0][1]), round(F1[i],2), round(Precision[i],2), round(Recall[i],2) ]
data.insert(i,x)
# write label and data
with open('Sup_Output.csv', 'w') as csvfile3:
datawriter = csv.writer(csvfile3, delimiter=',')
datawriter.writerow(label)
for data in data:
datawriter.writerow(data)