This notebook contains the code necessary to evaluate the performance of a convolutional neural network used in the false positive reduction step of a CAD system for pulmonary nodule detection. The set of candidates is given.
Two configurations are considered:
The two networks are evaluated on 1/5 of the publicly available LIDC-IDRI data set in a cross-validation fashion (although only one of the five folds is considered in this notebook).
The code in this notebook is based on the Theano library.
First, we import the libraries that we will be using in this example. Most of thse libraries can be installed by simply installing a distribution of Anaconda
In [ ]:
import numpy as np
import pylab as plt
import csv
import cPickle
import theano
import math
from theano import tensor as T
from theano.tensor.nnet.conv import conv2d
from theano.tensor.signal.downsample import max_pool_2d
from theano.sandbox.rng_mrg import MRG_RandomStreams as RandomStreams
from sklearn import metrics
%matplotlib inline
projectDir = '/home/user/nfbia15-deep-learning/examples/lungnodules/' # <-- set this local folder
srng = RandomStreams()
In [ ]:
def dropout(X, p=0.):
if p > 0:
retain_prob = 1 - p
X *= srng.binomial(X.shape, p=retain_prob, dtype=theano.config.floatX)
X /= retain_prob
return X
def softmax(X):
e_x = T.exp(X - X.max(axis=1).dimshuffle(0, 'x')) # 0 < e_x <= 1
return e_x / e_x.sum(axis=1).dimshuffle(0, 'x')
def rectify(X):
# ReLU non linearity
return T.maximum(X, 0.)
def init_weights(shape):
return theano.shared(floatX(np.random.randn(*shape) * 0.01))
def floatX(X):
return np.asarray(X, dtype=theano.config.floatX)
The function getDataset loads the data set from a numpy npz file. The file contains a field X for data and a field y for labels. Additionally, for the specific case of lung nodules, the position (in world coordinates) of the nodule and the seriesUID of the image can be loaded by setting the flag allInfo.
In [ ]:
def getDataset(datasetPath, allInfo = False):
npzFile = np.load(datasetPath)
X = npzFile['X']
y = npzFile['y']
if allInfo is True:
pos = npzFile['pos']
image = npzFile['img']
return X,y,pos,image
else:
return X,y
The Receiver Operator Characteristic curve (ROC) is a useful tool for evaluating the performance of a binary classification problem such as nodule detection (classification of nodules versus non nodules). In this notebook, we use the function implemented in the scikit-learn library. The input file is a csv file where the classification results are stored.
In [ ]:
def computeROC(inputfilename, figuretitle):
h = open(inputfilename,'rb')
csvreader = csv.reader(h)
ytest = []; preds = []
first = True
for row in csvreader:
if first:
groundtruth_idx = row.index('o_0')
prediction_idx = row.index('o_2')
first = False
else:
ytest.append(float(row[groundtruth_idx]))
preds.append(float(row[prediction_idx]))
# compute ROC curve
fpr, tpr, _ = metrics.roc_curve(ytest, preds)
auc = metrics.auc(fpr,tpr)
# figure
fig = plt.figure('ROC')
plt.title(figuretitle)
plt.xlim([0.0,1.0])
plt.ylim([0.0,1.0])
plt.gca().set_aspect(1.0)
plt.plot([0.0, 1.0],[0.0, 1.0],'k')
plt.xlabel('1 - specificity')
plt.ylabel('sensitivity')
ax = fig.gca(); ax.set_aspect('equal')
ax = fig.gca()
ax.set_xticks(np.arange(0.0,1.1,0.1))
ax.set_yticks(np.arange(0.0,1.1,0.1))
plt.grid()
plt.axis((0,1,0,1))
rocCurve, = plt.plot(fpr, tpr, linewidth=3)
plt.legend([rocCurve],['AUC = '+str(round(auc,4))[:5]+''],loc=4,prop={'size':12})
print 'AUC = '+str(auc)
The Free-Response ROC curve (FROC) shows the performance in terms of sensitivity as a function of the number of false positive per scan. It is a useful tool to understand the behaviour of the CAD system, as well as to set an operating point (by picking the threshold that corresponds to the desired amount of false positives per scan). The input for this function is the same csv file used for the ROC curve.
In [ ]:
def computeFROC(inputfilename, figuretitle, candidateDetectionSensitivity, datafilename, npoints=np.inf):
# collect number of cases
fh = open(inputfilename,'rb')
csvreader = csv.reader(fh)
# get indexes from the header
header = csvreader.next()
fh.seek(0)
seriesuid_idx = header.index('i_seriesuid')
groundtruth_idx = header.index('o_0')
prediction_idx = header.index('o_2')
seriesuids = []
thresholds = []
first = True
for csvrow in csvreader:
if first:
first = False
else:
seriesuids.append(csvrow[seriesuid_idx])
thresholds.append(float(csvrow[prediction_idx]))
fh.seek(0)
# remove duplicates
unique_seriesuids = list(set(seriesuids))
thresholds = sorted(list(set(thresholds)))
print str(len(unique_seriesuids))+' scans found'
# take npoints values of threholds
n = min(npoints, len(thresholds))
pts = range(0,len(thresholds),int(float(len(thresholds))/float(n)))
thresholds = [thresholds[i] for i in pts]
thresholds.append(1.0)
# fill in structured data
nlines = len(list(csvreader))-1
groundtruth_array = np.zeros((nlines, 1))
predictions_array = np.zeros((nlines, 1))
seriesuids_array = np.zeros((nlines, 1))
first = True
idx = 0
fh.seek(0)
for csvrow in csvreader:
if first:
first = False
else:
groundtruth_array[idx] = float(csvrow[groundtruth_idx])
predictions_array[idx] = float(csvrow[prediction_idx])
idx += 1
# assign an index to seriesuids
for idx in range(len(seriesuids)):
seriesuid = seriesuids[idx]
seriesuids_array[idx] = unique_seriesuids.index(seriesuid)
# compute actual values
nthresholds = len(thresholds)
sensitivities = np.zeros((nthresholds, 1))
fps_per_scan = np.zeros((nthresholds, 1))
nseriesuids = len(unique_seriesuids)
# compute froc
for th in range(nthresholds):
# set the threshold
threshold = thresholds[th]
positives = (predictions_array >= threshold)
# compute false positives and sensitivity
npositive = sum(groundtruth_array == 1.0)
tps_this_threshold = sum(positives[groundtruth_array == 1.0])
fps_this_threshold = sum(positives[groundtruth_array == 0.0])
sens_this_threshold = float(tps_this_threshold)/float(npositive)
sensitivities[th] = sens_this_threshold
fps_per_scan[th] = float(fps_this_threshold)/float(nseriesuids)
fps_per_scan = np.vstack((np.inf, fps_per_scan, -np.inf))
sensitivities = np.vstack((sensitivities[0], sensitivities, sensitivities[-1]))
# save results as csv file, used to compare the two curves later
ho = open(datafilename, 'w')
for i in range(len(fps_per_scan)):
ho.write(str(fps_per_scan[i][0])+','+str(candidateDetectionSensitivity*sensitivities[i][0])+'\n')
ho.close()
# plot curve
plt.figure()
plt.plot(fps_per_scan, candidateDetectionSensitivity*sensitivities, linewidth=3)
plt.xscale('log')
plt.ylim(0, 1)
plt.xlabel('Average number of false positives per scan')
plt.ylabel('Sensitivity')
plt.title('FROC Analysis')
plt.grid(b=True, which='both')
plt.tight_layout()
In this section, we will evaluate the performance of the convolutional network trained using only one 2D view of the noudule (the axial view). The network input is a 2D view of the nodule candidate, obtained by cropping the slice in the axial view with a bounding box of 50 mm x 50 mm. Patches are resized to 64x64 pixel. In this example, we will also define functions and variables necessary to train the network, but we will not train the network.
First, we point to the data on disk:
In [ ]:
modelsDir = projectDir+'models/' # where the pre-trained network is located
modelPath = modelsDir + '1View.pkl' # pre-trained network, it contains the parameters
outputPath = modelsDir + 'output_1_view.csv' # file containing classification results
print "model file: " + modelPath
We define the following network:
In [ ]:
# uses the variables to define the model
def model_1view(X, w1, w2, w3, w4, w_o, p_drop_hidden):
'''
'''
l1a = rectify(conv2d(X, w1, border_mode='valid')) #border_mode='full'
l1 = max_pool_2d(l1a, (2, 2))
l2a = rectify(conv2d(l1, w2))
l2 = max_pool_2d(l2a, (2, 2))
l3a = rectify(conv2d(l2, w3))
l3b = max_pool_2d(l3a, (2, 2))
l3 = T.flatten(l3b, outdim=2)
l4 = rectify(T.dot(l3, w4))
l4 = dropout(l4, p_drop_hidden)
pyx = softmax(T.dot(l4, w_o))
return l1, l2, l3, pyx, l1a, l2a, l3a
Before training, we need to initialize the parameters of the network. In this example, we only use weights without biases. The weights are initialized with random values by the function init_weights() defined above.
We also need to initialize a variable X which will represent our test data. The type ftensor4() is specific for the Theano library.
In [ ]:
# define weights and input data symbolic variables for 1 view convnet
w1_1view = init_weights((24, 1, 5, 5))
w2_1view = init_weights((32, 24, 3, 3))
w3_1view = init_weights((48, 32, 3, 3))
w4_1view = init_weights((48 * 6 * 6, 16))
w_o_1view = init_weights((16, 2))
X = T.ftensor4()
Now we pass the defined variable to the model function, and define the theano function that takes X as input and return the prediction, as well as the internal feature maps l1a, l2a, l3a and the weights w1, w2, w3, which we will use for visualization purpose later.
In [ ]:
l1_1view, l2_1view, l3_1view, py_x_1view, l1a_1view, l2a_1view, l3a_1view = model_1view(X, w1_1view, w2_1view, w3_1view, w4_1view, w_o_1view, 0.0)
y_x_soft_1view = py_x_1view # little abuse of notation
predictSoft_1view = theano.function(inputs=[X], outputs=[y_x_soft_1view, l1a_1view, l2a_1view, l3a_1view, w1_1view, w2_1view, w3_1view], allow_input_downcast=True)
In this example, we will not train the network, it would take too long. Since we already trained the network beforehand, we can load the trained weights and use the network to classify candidates. The following code is for loading the weights into memory. The pre-trained model has been saved as pkl file and can be loaded from disk.
In [ ]:
def loadModel_1view(modelPath):
saved_file = open(modelPath, 'rb')
w1_1view.set_value(cPickle.load(saved_file), borrow=True)
w2_1view.set_value(cPickle.load(saved_file), borrow=True)
w3_1view.set_value(cPickle.load(saved_file), borrow=True)
w4_1view.set_value(cPickle.load(saved_file), borrow=True)
w_o_1view.set_value(cPickle.load(saved_file), borrow=True)
In [ ]:
# load model
loadModel_1view(modelPath)
Once the weights are loaded, we can visualize them. For this purpose, we first define a function arrangeMaps, that takes a 4D matrix and re-arrange its content into a big square matrix, for better visualization.
In [ ]:
def arrangeMaps(x,nPerRow):
nchannels,nfilters,nr,nc = np.shape(x)
nmatrices = nchannels*nfilters
x = np.reshape(x,(nmatrices,nr,nc))
nrMatrix = int(nr*np.ceil((float(nchannels*nfilters)/float(nPerRow))))
ncMatrix = int(nPerRow*nc)
nPerCol = int(nrMatrix/nr)
xmap = np.inf*np.ones((nrMatrix + nPerCol-1, ncMatrix + nPerRow-1))
xmapIdx = 0
for i in range(0,nPerCol):
for j in range(0,nPerRow):
if xmapIdx < nmatrices:
xmap[i*nr+i:(i+1)*nr+i, j*nc+j:(j+1)*nc+j] = np.squeeze(x[xmapIdx,:,:])
xmapIdx += 1
return xmap
In [ ]:
plt.rcParams['figure.figsize'] = 6, 6
w1 = w1_1view.get_value().swapaxes(0,1)
print 'filter size: '+str(w1.shape)
print 'number of filters: '+str(w1.shape[0] * w1.shape[1])
w1_fig = arrangeMaps(w1,nPerRow=int(np.ceil(math.sqrt(w1.shape[0]*w1.shape[1]))))
plt.imshow(w1_fig, cmap='gray', interpolation='nearest')
plt.axis('off')
Values in white are positive, values in black are negative, values in gray are (close to) zero.
In [ ]:
plt.rcParams['figure.figsize'] = 16, 12
w2 = w2_1view.get_value().swapaxes(0,1)
print 'filter size: '+str(w2.shape)
print 'number of filters: '+str(w2.shape[0] * w2.shape[1])
w2_fig = arrangeMaps(w2,nPerRow=int(np.ceil(math.sqrt(w2.shape[0]*w2.shape[1]))))
plt.imshow(w2_fig, cmap='gray', interpolation='nearest')
plt.axis('off')
In [ ]:
plt.rcParams['figure.figsize'] = 16, 12
w3 = w3_1view.get_value().swapaxes(0,1)
print 'filter size: '+str(w3.shape)
print 'number of filters: '+str(w3.shape[0] * w3.shape[1])
w3_fig = arrangeMaps(w3,nPerRow=int(np.ceil(math.sqrt(w3.shape[0]*w3.shape[1]))))
plt.imshow(w3_fig, cmap='gray', interpolation='nearest')
plt.axis('off')
Now we can use the pre-trained network to classify the data set containing one view for each candidate. For this example, we used ~1/5 of the LIDC-IDRI data set, in a five-fold cross-validation fashion. The pre-trained network was trained using candidated in the other 4 folds, 3 folds for training, 1 fold for validation, and some data augmentation was used.
The detection of candidates is part of an ongoing research project, which gives a candidate detection sensitivity of 0.944. In this example, we will not run rhe network on the whole test data set (>2GB), since it takes several minutes. We did it beforehand, and we stored the classification output in a csv file (output_1_view.csv), that we can load to compute the ROC and FROC curves.
The function classifyDataset_1plane takes a dataset file and a pre-trained network as input, and classifies the candidates. The results are stored as a csv file.
In [ ]:
def classifyDataset_1plane(datasetPath, modelPath, outputPath, showFig=False, nExamples=np.inf):
dsX,dsY,dsPos,dsImage = getDataset(datasetPath, allInfo = True)
dsX = dsX.reshape(-1, 1, 64, 64)
res = open(outputPath,'w')
res.write('f_0,o_0,o_1,o_2,i_seriesuid,i_worldVectorX,i_worldVectorY,i_worldVectorZ\n')
for i in range(min(len(dsX), nExamples)):
if i % 1000 == 0:
print "Progress = " + str(100*float(i)/len(dsX)) + "%"
xcoord = dsPos[i][0]
ycoord = dsPos[i][1]
zcoord = dsPos[i][2]
img = dsImage[i]
y = dsY[i][1]
imarray = [dsX[i,:,:,:]]
pred,l1a,l2a,l3a,w1,w2,w3 = predictSoft_1view(imarray)
# show image with classification steps
if showFig:
showClassification(imarray,[l1a,l2a,l3a],[w1,w2,w3],pred,y)
res.write('0,'+str(dsY[i][1])+','+str(pred[0][0])+','+str(pred[0][1])+','+img+','+str(xcoord)+','+str(ycoord)+','+str(zcoord)+'\n')
res.close()
To better uderstand the performance of the network, we visualize the input patch, the three internal feature maps and the final prediction. For this purpose, we made a small data set containing 100 nodules and 100 false positives. By the default, the number of examples to show is nExamples=10, but you can increase it up to nExamples=200.
In [ ]:
def showClassification(imarray,featMaps,weights,prediction,y):
nFeatMaps = len(featMaps)
img = np.squeeze(imarray)
# show input image
plt.subplot(1,nFeatMaps+2,1)
plt.imshow(img, cmap='gray')
plt.axis('off')
if y == 0:
plt.title('non nodule')
else:
plt.title('nodule')
# show feature maps
for i in range(nFeatMaps):
plt.subplot(1,nFeatMaps+2,2+i)
fmap = arrangeMaps(featMaps[i],nPerRow=int(np.ceil(math.sqrt(featMaps[i].shape[0] * featMaps[i].shape[1]))))
plt.imshow(fmap, cmap='jet',interpolation='nearest')
plt.axis('off')
plt.title('C'+str(i+1)+'')
# show prediction
plt.subplot(1,nFeatMaps+2,nFeatMaps+2)
pos = range(len(prediction.flatten()))
plt.bar(pos, prediction.flatten(), width=1)
plt.xticks([0.5,1.5], ('non nodule', 'nodule'))
plt.title('prediction')
plt.show()
# classify
datasetDir = projectDir+'/data/'
modelsDir = projectDir+'/models/'
testPath = datasetDir + '1ViewTestset_small.npz'
modelPath = modelsDir + '1View.pkl'
outputPath = modelsDir + 'output_1_view_small.csv'
plt.rcParams['figure.figsize'] = 16, 2.5
classifyDataset_1plane(testPath, modelPath, outputPath, showFig=True, nExamples=10)
In [ ]:
# compute ROC curve
plt.rcParams['figure.figsize'] = 6, 6
inputfile = modelsDir + 'output_1_view.csv'
computeROC(inputfile, 'ROC - one patch per nodule')
Computing the FROC curve takes a bit more than the ROC curve. To speed up the computation, you can reduce the value in npoints, which indicates how many threshold values we consider to compute FPs/scan.
Note that the value of 'candidateDetectionSensitivity' is also passed as a parameter. We have to take into account for the candidate detection sensitivity, which is the maximum we can get, because some nodules in the scan might not be in the set of candidates. If we did not consider it, the FROC curve will get up to a sensitivity of 1.0, since all the nodule candidates will be detected.
In [ ]:
# compute FROC
plt.rcParams['figure.figsize'] = 8, 5
inputfile = modelsDir + 'output_1_view.csv'
figurefilename = projectDir+'froc_1_view.pdf'
figuretitle = 'FROC - one patch per nodule'
candidateDetectionSensitivity = 0.944 # pre-computed
datafilename = projectDir+'froc_data_1_view.csv'
computeFROC(inputfile, figuretitle, candidateDetectionSensitivity, datafilename, npoints=100)
We repeat the same procedure but for a network that receives three views of the nodule simultaneously, in particular, three patches along the axial, coronal and sagittal view are considered.
The architecture consists of three streams of convolutional and max-pooling layers, which share the weights. The final softmax layer is connected to the three fully-connected layers and makes the final prediction.
In [ ]:
modelsDir = projectDir+'models/'
modelPath = modelsDir + '3Views.pkl'
outputPath = modelsDir + 'output_3_views.csv'
rocPath = modelsDir + 'roc_3_views.csv'
print "model file: " + modelPath
The weight configuration is exactly the same as the network working on 1 view.
In [ ]:
# define weights and input data symbolic variables for 3 views convnet
w1_3views = init_weights((24, 1, 5, 5))
w2_3views = init_weights((32, 24, 3, 3))
w3_3views = init_weights((48, 32, 3, 3))
w4_3views = init_weights((48 * 6 * 6, 16))
w_o_3views = init_weights((16, 2))
X1 = T.ftensor4()
X2 = T.ftensor4()
X3 = T.ftensor4()
Here we define the model to work on the three views. The weights are shared, so we use w1 for the first layer in the three streams, w2 for the second layer in the three streams etc.
In [ ]:
def model_3views(X1, X2, X3, w1, w2, w3, w4, w_o, p_drop_hidden):
# layer 1
l1a_p1 = rectify(conv2d(X1, w1, border_mode='valid'))
l1_p1 = max_pool_2d(l1a_p1, (2, 2))
l1a_p2 = rectify(conv2d(X2, w1, border_mode='valid'))
l1_p2 = max_pool_2d(l1a_p2, (2, 2))
l1a_p3 = rectify(conv2d(X3, w1, border_mode='valid'))
l1_p3 = max_pool_2d(l1a_p3, (2, 2))
# layer 2
l2a_p1 = rectify(conv2d(l1_p1, w2))
l2_p1 = max_pool_2d(l2a_p1, (2, 2))
l2a_p2 = rectify(conv2d(l1_p2, w2))
l2_p2 = max_pool_2d(l2a_p2, (2, 2))
l2a_p3 = rectify(conv2d(l1_p3, w2))
l2_p3 = max_pool_2d(l2a_p3, (2, 2))
# layer 3
l3a_p1 = rectify(conv2d(l2_p1, w3))
l3b_p1 = max_pool_2d(l3a_p1, (2, 2))
l3_p1 = T.flatten(l3b_p1, outdim=2)
l3a_p2 = rectify(conv2d(l2_p2, w3))
l3b_p2 = max_pool_2d(l3a_p2, (2, 2))
l3_p2 = T.flatten(l3b_p2, outdim=2)
l3a_p3 = rectify(conv2d(l2_p3, w3))
l3b_p3 = max_pool_2d(l3a_p3, (2, 2))
l3_p3 = T.flatten(l3b_p3, outdim=2)
# layer 4
l4_p1 = rectify(T.dot(l3_p1, w4))
l4_p2 = rectify(T.dot(l3_p2, w4))
l4_p3 = rectify(T.dot(l3_p3, w4))
# concatenate the outputs of 3 fully connected layers
l4 = T.concatenate([l4_p1, l4_p2, l4_p3], axis=1)
l4 = dropout(l4, p_drop_hidden)
pyx = softmax(T.dot(l4, w_o))
return l1_p1, l2_p1, l3_p1, pyx, l1a_p1, l2a_p1, l3a_p1
In [ ]:
# 3 views
l1_p1_3views, l2_p1_3views, l3_p1_3views, py_x_3views, l1a_3views, l2a_3views, l3a_3views = model_3views(X1, X2, X3, w1_3views, w2_3views, w3_3views, w4_3views, w_o_3views, 0.0)
y_x_soft_3views = py_x_3views
predictSoft_3views = theano.function(inputs=[X1, X2, X3], outputs=[y_x_soft_3views, l1a_3views, l2a_3views, l3a_3views, w1_3views, w2_3views, w3_3views], allow_input_downcast=True)
In [ ]:
def loadModel_3views(modelPath):
save_file = open(modelPath, 'rb')
w1_3views.set_value(cPickle.load(save_file), borrow=True)
w2_3views.set_value(cPickle.load(save_file), borrow=True)
w3_3views.set_value(cPickle.load(save_file), borrow=True)
w4_3views.set_value(cPickle.load(save_file), borrow=True)
w_o_3views.set_value(cPickle.load(save_file), borrow=True)
In [ ]:
loadModel_3views(modelPath)
In [ ]:
def classifyDataset_3planes(datasetPath, modelPath, outputPath):
dsX,dsY,dsPos,dsImage = getDataset(datasetPath, allInfo = True)
dsX1 = dsX[:,:,:,0]
dsX1 = dsX1.reshape(-1, 1, 64, 64)
dsX2 = dsX[:,:,:,1]
dsX2 = dsX2.reshape(-1, 1, 64, 64)
dsX3 = dsX[:,:,:,2]
dsX3 = dsX3.reshape(-1, 1, 64, 64)
res = open(outputPath,'w')
res.write('f_0,o_0,o_1,o_2,i_seriesuid,i_worldVectorX,i_worldVectorY,i_worldVectorZ\n')
for i in range(len(dsX)):
if i % 100 == 0:
print "Progress = " + str(100*float(i)/len(dsX)) + "%"
xcoord = dsPos[i][0]
ycoord = dsPos[i][1]
zcoord = dsPos[i][2]
img = dsImage[i]
imarray1 = [dsX1[i,:,:,:]]
imarray2 = [dsX2[i,:,:,:]]
imarray3 = [dsX3[i,:,:,:]]
pred = predictSoft_3views(imarray1,imarray2,imarray3)[0]
res.write('0,'+str(dsY[i][1])+','+str(pred[0][0])+','+str(pred[0][1])+','+img+','+str(xcoord)+','+str(ycoord)+','+str(zcoord)+'\n')
res.close()
We do not classify the data set here, as for the example with 1 view, we did it beforehand and saved the classification results in a csv file. Now we evaluate the performance via ROC and FROC analysis.
In [ ]:
# compute ROC
plt.rcParams['figure.figsize'] = 6, 6
inputfile = outputPath
figurefilename = projectDir+'roc_3_views.pdf'
figuretitle = 'ROC - three patches per nodule'
computeROC(inputfile, figuretitle)
In [ ]:
# compute FROC
plt.rcParams['figure.figsize'] = 8, 5
figurefilename = projectDir+'froc_3_views.pdf'
figuretitle = 'FROC - three patches per nodule'
candidateDetectionSensitivity = 0.944 # pre-computed
datafilename = projectDir+'froc_data_3_views.csv'
computeFROC(inputfile, figuretitle, candidateDetectionSensitivity, datafilename, npoints=100)
In [ ]:
plt.rcParams['figure.figsize'] = 8, 5
datafilename1 = modelsDir+'froc_data_1_view.csv'
datafilename3 = modelsDir+'froc_data_3_views.csv'
froc1 = open(datafilename1,'r').readlines()
froc3 = open(datafilename3,'r').readlines()
x1 = []; y1 = []
x3 = []; y3 = []
for line in froc1:
xy = line.rstrip().split(',')
x1.append(float(xy[0])); y1.append(float(xy[1]))
for line in froc3:
xy = line.rstrip().split(',')
x3.append(float(xy[0])); y3.append(float(xy[1]))
plt.figure()
froc1, = plt.plot(x1, y1, linewidth=3, color='red')
froc3, = plt.plot(x3, y3, linewidth=3, color='green')
plt.xscale('log')
plt.yscale('linear')
plt.xlabel('Average number of false positives per scan')
plt.ylabel('Sensitivity')
plt.title('FROC Analysis')
plt.grid(b=True, which='both')
plt.tight_layout()
plt.legend([froc1, froc3],['1 view','3 views'],loc=4,prop={'size':12})
plt.show()