In [52]:
#requires from python: matplotlib, numpy and scipy
# from debian: ffmpeg and libx264
import autoPlanLib as myTools #this has some custom functions to clean up the code here
In [74]:
#load patient DVH data (settings are in autoPlanLib.py)
patients = []
# read ting's data
# for PTV DVH line pairs are dose and volume (for 24 patients)
f = open('ting_data/DVH_PTV_23Patients.txt')
parity = 0 #tracks line pairs
for line in f.readlines():
if len(line) >1:
if parity == 0:
#we have a new pair of dose,volume lines (and a new patient)
patients.append([])
patients[-1].append(myTools.Organ())
#start with dose
patients[-1][-1].name = 'PTV'
patients[-1][-1].dose = map(lambda x: float(x), line.rstrip().split(' '))
patients[-1][-1].n = len(patients[-1][-1].dose)
#increment parity
parity = 1
else:
#this is volume data for the above dose data
patients[-1][-1].volume = map(lambda x: float(x)/100., line.rstrip().split(' '))
#increment parity
parity = 0
#now read in the Rectum-DWH
f = open('ting_data/Rectum-DWHData_23Patients.txt')
#the first line is doses:
doses = map(lambda x: float(x), f.readline().rstrip().split(' '))
pid = 0
for line in f.readlines():
if len(line.rstrip()) >1:
#this is wall coverage
patients[pid].append(myTools.Organ())
patients[pid][-1].name = 'Rectum DWH'
patients[pid][-1].n = len(doses)
patients[pid][-1].dose = doses
patients[pid][-1].volume = map(lambda x: float(x)/100., line.rstrip().split(' '))
pid += 1
f.close()
print "total patients:%i" % pid
f.close()
print "total patients:%i" % len(patients)
#now read in the bladder-DWH
f = open('ting_data/Bladder-DWHData_23Patients.txt')
#the first line is doses:
doses = map(lambda x: float(x), f.readline().rstrip().split(' '))
pid = 0
for line in f.readlines():
if len(line.rstrip()) >1:
#this is wall coverage
patients[pid].append(myTools.Organ())
patients[pid][-1].name = 'Bladder DWH'
patients[pid][-1].n = len(doses)
patients[pid][-1].dose = doses
patients[pid][-1].volume = map(lambda x: float(x)/100., line.rstrip().split(' '))
pid += 1
f.close()
print "total patients:%i" % pid
In [75]:
myTools = reload(myTools) #when needed
In [76]:
#Some useful parameters
Rx_dose = 100
num_pat = size(patients,0) #number of patients in library (60)
num_orgs = 3 #number of organs #len(patients[pat])
In [77]:
#plotting DVH of library
num_pat = 23
fig, ax = plt.subplots(1,3,figsize=(18,5))
for pat in range(num_pat):
for org in range(num_orgs):
data = patients[pat][org]
ax[org].plot(data.dose,multiply(data.volume,100.0))#,organ_fmt[o])
if pat == num_pat - 1:
if org == 0:
ax[org].set_title("All " + data.name)
ax[org].set_xlabel('Dose (cGy)')
ax[org].set_ylabel('Volume (%)')
ax[org].set_ylim(0,100)
else:
ax[org].set_title("All " + data.name)
ax[org].set_xlabel('Dose (cGy)')
ax[org].set_ylabel('Area (%)')
ax[org].set_ylim(0,100)
Here we have some freedom in how we re-sample the data for use in PCA. PCA requires a data matrix, each row being a observation (patient), and each column representing some measurement. The PCA library we use also assumes that the number of columns in the data matrix is less than or equal to the number or rows.
A few different sampling schemes were tested, but in the end they all gave very similar results.
Below is a scheme to weight the sampling grid quadratically towards lower volume (%) values (higher dose regions).
In [6]:
# say we want to weight samples by x^2
# integrate, then invert we get 3/x^3, then uniformly sample
n_pts = 20
grid_points = zeros((num_orgs,n_pts))
g = mgrid[0.0:1.0:(1./float(n_pts))] #a uniform grid with desired number of sampling points
sampling = map(lambda x: pow(x,1.0/3.0)*100, g) #map uniform grid through weighting function
sampling = abs(array(sampling)-98) #here the sampling points are effectively reversed and shifted
print "Volume sampling points (%):"
print sampling
#Vols for PTV, RECT and BLAD:
for i in range(3):
grid_points[i] = sampling / 100.0
Below are 2 hand picked sampling scheme which treats the PTV and OARS differently
In [7]:
n_pts = 12 #points sampled from each organ
grid_points = zeros((num_orgs,n_pts))
#Vols for PTV:
grid_points[0] = array( (99,98,97,96,95,90,50, 5, 4, 3, 2, 1) ) / 100.0
#Vols for RECT:
grid_points[1] = array( (90,70,60,50,40,30,20,10, 9, 7, 5, 1) ) / 100.0
#Vols for BLAD:
grid_points[2] = array( (90,70,60,50,40,30,20,10, 9, 7, 5, 1) ) / 100.0
print "Volume sampling points (%):"
print 100*grid_points
In [23]:
#another hand picked sampling scheme with fiewer points (for simplicity)
n_pts = 6 #points sampled from each organ
grid_points = zeros((num_orgs,n_pts))
#Vols for PTV:
grid_points[0] = array( (99,95,90,70,50,2) ) / 100.0
#Vols for RECT:
grid_points[1] = array( (90,60,30,20,10,5) ) / 100.0
#Vols for BLAD:
grid_points[2] = array( (90,60,30,20,10,5) ) / 100.0
print "Volume sampling points (%):"
print 100*grid_points
In [56]:
#a uniform sampling scheme with 20 sampling points per organ
n_pts = 7
low_vol = 1.
high_vol = 99.
a = frange(low_vol, high_vol, (high_vol-low_vol)/float(n_pts-1))
a = a[::-1] #reversed
grid_points = repeat([a/100.],3,axis=0)
In [57]:
#Plot the grid
grid_points100 = 100*grid_points #this gives %dose units
plot(grid_points100[0],grid_points100[0],'o')
title('Sampling points for Volume (%)')
show()
In [78]:
gridDim = len(grid_points[0])
V_grid = zeros((num_pat,num_orgs*gridDim))
for pat in range(num_pat):
for org in range(num_orgs):
data = patients[pat][org]
for i in range(len(grid_points[org])):
V_grid[pat,org*gridDim+i]=myTools.DoseAtVol(data.dose,data.volume,grid_points[org][i])/Rx_dose
#this defines the indices for each organ in the rows of V_grid
r = (range(0,n_pts),range(n_pts,2*n_pts),range(2*n_pts,3*n_pts))
In [79]:
fig, axes = plt.subplots(1,3,figsize=(18,4))
for pat in range(num_pat):
axes[0].plot( V_grid[pat,r[0]], grid_points100[0],'-m', alpha=.3)
axes[1].plot( V_grid[pat,r[1]], grid_points100[1],'-r', alpha=.3)
axes[2].plot( V_grid[pat,r[2]], grid_points100[2],'-b', alpha=.3)
#axes[i,2].plot([],[],'--k',label='Neg(%i)'% m_idx)
#axes.plot(V_grid[p,:],'.',alpha = 0.5)
axes[0].set_title(patients[0][0].name)
map(lambda x: x.set_xlabel('Dose/RxDose'),axes)
axes[0].set_ylabel('Volume (%)')
axes[1].set_title(patients[0][1].name)
axes[2].set_title(patients[0][2].name)
show()
In [80]:
from matplotlib.mlab import PCA
V_pca = PCA(V_grid)
scl = 1.2;
figure(figsize=(scl*7,scl*5))
title('PCA Fractions')
plot(V_pca.fracs*100,'or') # the proportion of variance of each of the principal components
ylabel("Variance (%)")
xlabel("Principal Component")
ylim(-1,46)
xlim(-1.0,20.5)
show()
In [81]:
#scores for the first 3 modes for each patient
print"PC0 PC1 PC2 PC3 PC4 PC5"
for p in range(len(patients)):
print "%f %f %f %f %f %f" % (V_pca.Y[p,0],V_pca.Y[p,1],V_pca.Y[p,2],V_pca.Y[p,3],V_pca.Y[p,4],V_pca.Y[p,5])
The plot below gives a hint as to how the data matrix is sorted. Each row represents a patient DVH and the first 1/3 of the columns sample the PTV, the middle 1/3 sample the Rectum, and the last 1/3 sample the Bladder.
The PCA modes represent (unit) directions in the high dimensional data space (after the mean for each column has been subtracted and the standard deviation is divided out).
In [82]:
n_comps = 6
p = V_pca
fig, axes = plt.subplots(figsize=(18,7))
for i in range(n_comps):
axes.plot(p.Wt[i,:],label='PC %i'%i)
axes.plot(zeros(num_orgs*n_pts),':k')
axes.set_xlabel('Volume (%) for PTV (Purple), Rectum (Red), Bladder (Blue)')
axes.set_ylabel('weights/correlations (arb)')
axes.set_title('First %i PCA modes' % n_comps)
axes.legend(loc=0)
#format Ticks
axes.set_xticks(range(3*n_pts))
#make integer axis of Vol axis
l = axes.set_xticklabels(array(concatenate((grid_points100[0],grid_points100[1],grid_points100[2])),int)) #output hidden
junk = map(lambda x: x.set_color('purple'), l[:n_pts] )
junk = map(lambda x: x.set_color('red'), l[n_pts:2*n_pts] )
junk = map(lambda x: x.set_color('blue'), l[2*n_pts:] )
In [83]:
n_comps = 6
p = V_pca
fig2, axes2 = plt.subplots(n_comps,3,figsize=(18,n_comps*6))
for i in range(n_comps):
axes2[i,0].plot( zeros((len(grid_points[0]),1)),grid_points100[0],':k')
axes2[i,0].plot(p.Wt[i,r[0]],grid_points100[0],'-m',label="PTV")
axes2[i,0].plot(p.Wt[i,r[1]],grid_points100[1],'-r',label="Rectum")
axes2[i,0].plot(p.Wt[i,r[2]],grid_points100[2],'-b',label="Bladder")
axes2[i,0].set_title('Principal Componant Vector # %i' % i)
axes2[i,0].set_ylabel('Volume (%)')
axes2[i,0].set_xlabel('Weight/Correlation (arb)')
axes2[i,0].legend(loc=0)
#plot PC dimension (componants for each patient)
axes2[i,1].plot( p.Y[:,i],p.Y[:,i],'sr',mew=0,ms=20,alpha = 0.15)
axes2[i,1].set_title('PC %i vs. PC %i (squares are patients)' % (i,i))
axes2[i,1].set_xlabel('Weight/Correlation (arb)')
#Find Extremum Patients for PC mode (Most Negative, Closest to Zero, and Most Positive)
neg_idx = find(min(p.Y[:,i]) == p.Y[:,i])
tempArray = map(lambda x: abs(x), p.Y[:,i]) #map absolute value
zero_idx = find(min(tempArray) == tempArray)
pos_idx = find(max(p.Y[:,i]) == p.Y[:,i])
c_fmt = ('m','r','b') #line colors for PTV, Rect, and Blad
axes2[i,2].plot([],[],'--k',label='Neg(%i)'% neg_idx)
axes2[i,2].plot([],[],':k',label='Zero(%i)' % zero_idx)
axes2[i,2].plot([],[],'-k',label='Pos(%i)' % pos_idx)
for org in range(3):
axes2[i,2].plot( patients[neg_idx][org].dose, 100.*array(patients[neg_idx][org].volume),'--'+c_fmt[org])
axes2[i,2].plot( patients[zero_idx][org].dose, 100.*array(patients[zero_idx][org].volume),':'+c_fmt[org])
axes2[i,2].plot( patients[pos_idx][org].dose, 100.*array(patients[pos_idx][org].volume),'-'+c_fmt[org])
axes2[i,2].set_xlim((0,Rx_dose*1.2))
axes2[i,2].set_ylim((0,110))
axes2[i,2].set_title('Extremum DVH along PC %i ' % i)
axes2[i,2].set_xlabel('Dose (cGy)')
axes2[i,2].set_ylabel('Volume (%)')
axes2[i,2].legend(loc='upper center',ncol=3,frameon=False)
In [84]:
#3D scatter in PCV (0,1,3) space
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(1, 1, 1, projection='3d')
pY = V_pca.Y
p = ax.scatter(pY[:,0], pY[:,1], pY[:,2],s=60,c='r',marker='o',edgecolors='none',alpha = 0.5)
#maybe try animation later
In [85]:
p = V_pca
proj = p.project(V_grid[0,:])
temp = zeros(proj.shape)
errors = zeros(proj.shape)
n_pats = len(patients) #number of patients to consider
n_pc = 6 #number of pca modes to include
#placeholders to find max error patient
max_error = 0.0
max_error_row = []
max_error_pca = []
for pat in range(n_pats):
row = V_grid[pat,:]
proj = p.project(row)
temp = zeros(proj.shape)
for i in range(n_pc):
temp = temp + multiply(p.Wt[i,:],proj[i])
#compute error between pca reduced representation
temp_err = abs(temp*p.sigma + p.mu - row)
#compute max error and compare to global
if max(temp_err) > max_error:
#save patient index
max_error_idx = pat
max_error_row = row.copy()
max_error_pca = temp.copy()
errors = errors + temp_err
In [86]:
#plot PCA errors
fig, ax = plt.subplots(1,2,figsize=(16,7))
ax[0].plot(100*errors[r[0]]/float(n_pats),grid_points100[0],':om',label="PTV")
ax[0].plot(100*errors[r[1]]/float(n_pats),grid_points100[0],':or',label="Rectum")
ax[0].plot(100*errors[r[2]]/float(n_pats),grid_points100[0],':ob',label="Bladder")
ax[0].set_title('Average Error with first %i PCA modes' % n_pc)
ax[0].legend(loc='best')
ax[0].set_ylabel('Volume (%)')
ax[0].set_xlabel('Absolute Error (% of RxDose)')
#example Origional and Reduced
dvh_pca = max_error_pca*p.sigma + p.mu
ax[1].plot(Rx_dose * dvh_pca[r[0]],grid_points100[0],'-m',label="PTV")
ax[1].plot(Rx_dose * dvh_pca[r[1]],grid_points100[1],'-r',label="Rectum")
ax[1].plot(Rx_dose * dvh_pca[r[2]],grid_points100[2],'-b',label="Bladder")
#plot(max_error_pca*p.sigma + p.mu)
ax[1].plot(Rx_dose * max_error_row[r[0]],grid_points100[0],'--m')
ax[1].plot(Rx_dose * max_error_row[r[1]],grid_points100[1],'--r')
ax[1].plot(Rx_dose * max_error_row[r[2]],grid_points100[2],'--b')
ax[1].plot([],[],'--k',label='Orig.')
ax[1].plot([],[],'-k',label='Proj.')
ax[1].set_ylabel('Volume (%)')
ax[1].set_xlabel('Dose (cGy)')
ax[1].legend(loc=0,ncol=2)
ax[1].set_title('Worst Max Error (Patient %i)' % max_error_idx)
show()