In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [206]:

    

from scipy.io import loadmat, savemat
from util import mosaic, qimshow, status_check
import matplotlib.pyplot as plt
from numpy.ma import masked_array
from skimage.filter import threshold_otsu
import time
import dcemri
import dicom


# SCRIPT FLAGS
plotting = True
reduce_problem = False
Rel = 4.5   # Relaxivity of Gd-DTPA at 3 T [s^-1 [mmol Gd-DTPA]^{-1}]   
flip_angle = 30 * pi / 180.0 # rad
TR = 5e-3   # sec
nx = 80
ny = 50
nt = 61
t = 6.0*arange(nt)  # 6 s per dynamic


data_dir = 'qiba/v9'
#file_ext = 'QIBA_v9_Tofts_S0_500_6s_0s_sigma_5'
file_ext = 'QIBA_v9_Tofts_S0_10000_6s_0s_sigma_100'


    

In [207]:

    

# add cubehelix3 color map
mpl.cm.register_cmap(name='cubehelix3', data=mpl._cm.cubehelix(gamma=1.0, s=0.4, r=-0.5, h=1.5))


    

In [208]:

    

data_dicom = zeros((nt, nx, ny))
for k in range(nt):
    file_name = '%s/%s_%03d.dcm' % (data_dir, file_ext, k+1)
    dcm = dicom.read_file(file_name)
    data_dicom[k,:,:] = dcm.pixel_array.astype('float')
    #print k, dcm.pixel_array.shape, file_name
    
data_dce = data_dicom[:,10:70,:]
nt, nx, ny = data_dce.shape
T1map = ones((nx, ny)) # s
#S0map = ones((nx, ny)) * 50000.0 #
data_aif = mean(mean(data_dicom[:,70:,:], axis=2), axis=1)


    

In [209]:

    

imshow(data_dce[:,:,:].max(axis=0), interpolation='nearest')


    

    
Out[209]:





<matplotlib.image.AxesImage at 0x1110a14d0>






    

In [210]:

    

plot(data_dce[:,-1,-1],'k.-')
title('example tissue signal curve')


    

    
Out[210]:





<matplotlib.text.Text at 0x1128c2c90>






    

In [211]:

    

# turn Sb into Cp
T1p = 1.440
R1p = 1.0 / T1p
Hct = 0.45
S0 = data_aif[:10].mean()
R1_eff_aif = dcemri.dce_to_r1eff(data_aif, S0, R1p, TR, flip_angle)
Cb = dcemri.r1eff_to_conc(R1_eff_aif.flatten(), R1p, Rel)  
Cp = Cb.flatten() / (1.0 - Hct)


    

    




converting DCE signal to effective R1
converting effective R1 to tracer tissue concentration

In [212]:

    

if plotting:
    figure(1)
    clf()
    plot(t, Cp, 'k.-')
    xlabel('time (s)')
    ylabel('$C_p$ (mM)')
    title('derived AIF')
    savefig('qiba/qiba_v9_aif.pdf')


    

In [213]:

    

data_dce = transpose(data_dce, (1,2,0)).copy('C')

if reduce_problem:
    nxx = 6
    nyy = 5
    data_dce_red = zeros((nxx, nyy, nt))
    t1_map_red = zeros((nxx, nyy))
    s0_map_red = zeros((nxx, nyy))
    for x in range(nxx):
        for y in range(nyy):
            tmp = data_dce[10*x:10*(x+1), 10*y:10*(y+1), :]
            data_dce_red[x,y,:] = mean(mean(tmp, axis=0), axis=0)
            tmp = T1map[10*x:10*(x+1), 10*y:10*(y+1)]
            t1_map_red[x,y] = mean(tmp)
    data_dce = data_dce_red
    T1map = t1_map_red
    nx = nxx
    ny = nyy


    

In [214]:

    

R1map = 1.0 / T1map


    

In [215]:

    

data_dce = reshape(data_dce, (nx, ny, nt))
S0 = data_dce[:,:,:10].mean(axis=2)
R1_eff = dcemri.dce_to_r1eff(data_dce, S0, R1map.T, TR, flip_angle)


    

    




converting DCE signal to effective R1

In [216]:

    

# show map of R1_eff 
R1_eff = reshape(R1_eff, (nx, ny, nt))
if plotting:
    figure(7)
    clf()
    imshow(mosaic(R1_eff.max(axis=2)), interpolation='nearest', cmap='cubehelix3')
    title('max $R_1(t) (s^{-1})$')
    xlabel('ro')
    ylabel('pe')
    colorbar()


    

In [217]:

    

# convert effecitve R1 to tissue concentration Ct
Ct = dcemri.r1eff_to_conc(transpose(R1_eff, (2,0,1)), R1map, Rel)
Ct = transpose(Ct, (1,2,0))


    

    




converting effective R1 to tracer tissue concentration

In [218]:

    

# show map of Ct
Ct = reshape(Ct, (nx, ny, nt))
if plotting:
    figure(9)
    clf()
    imshow(Ct.max(axis=2), interpolation='nearest', cmap='cubehelix3')
    title('max $C_t$ [mM Gd-DTPA]')
    colorbar()
    savefig('qiba/qiba_v9_Ct.pdf')
Ct = reshape(Ct, (-1, nt))


    

In [232]:

    

dcemri = reload(dcemri)
params, covs =  dcemri.fit_tofts_model(Ct, Cp, t, extended=False)


    

    




fitting perfusion parameters
using Standard Tofts-Kety
fitting 3000 voxels
10% complete, 1012 of 1125 s remain
20% complete, 938 of 1172 s remain
30% complete, 763 of 1090 s remain
40% complete, 631 of 1052 s remain
50% complete, 513 of 1027 s remain
60% complete, 393 of 983 s remain
70% complete, 287 of 959 s remain
80% complete, 182 of 910 s remain
90% complete, 87 of 870 s remain
100% complete, 0 of 837 s remain
837 s elapsed

In [233]:

    

Kt = zeros(nx*ny)
idxs_dce = range(nx*ny)
Kt[idxs_dce] = params[0]
ve = zeros(nx*ny)
ve[idxs_dce] = params[1]
Kt = reshape(Kt, (nx, ny))
ve = reshape(ve, (nx, ny))


    

In [234]:

    

# convert Kt to min$^{-1}$
Kt = Kt * 60.0


    

In [235]:

    

# PLOT RESULTS
if plotting:
    figure(13)
    clf()
    imshow(Kt, interpolation='nearest', cmap='cubehelix3', vmin=0, vmax=0.5)
    #plot(Kt.T,'o-')
    colorbar()
    title('$K^{trans}$ (min$^{-1}$)')
    savefig('qiba/qiba_v9_Kt.pdf')

    figure(14)
    clf()
    imshow(ve, interpolation='nearest', cmap='cubehelix3', vmin=0, vmax=1)
    #plot(ve,'o-')
    title('$v_e$')
    colorbar()
    savefig('qiba/qiba_v9_ve.pdf')


    

In [236]:

    

if reduce_problem:
    u = ones(1)
else:
    u = ones(10)
Kt_truth = hstack([0.01*u, 0.02*u, 0.05*u, 0.1*u, 0.2*u, 0.35*u])
ve_truth = hstack([0.01*u, 0.05*u, 0.1*u, 0.2*u, 0.5*u])
Kt_error = ((Kt.T - Kt_truth) / Kt.T).T
ve_error = (ve - ve_truth) / ve


    

In [237]:

    

if plotting:
    figure(10)
    imshow(Kt_error, interpolation='nearest', cmap='bwr', vmin=-2, vmax=2)
    title('relative error in $K^{trans}$')
    colorbar()
    savefig('qiba/qiba_v9_error_Kt.pdf')

    
    figure(11)
    imshow(ve_error, interpolation='nearest', cmap='bwr', vmin=-2, vmax=2)
    title('relative error in $v_e$')
    colorbar()
    savefig('qiba/qiba_v9_error_ve.pdf')


    

In [238]:

    

# mean relative error
print 'rms relative error in Kt:', abs(Kt_error).mean()
print 'rms relative error in ve:', abs(ve_error).mean()


    

    




rms relative error in Kt: 75551.4717183
rms relative error in ve: 372486.642768

In [239]:

    

h = hist(Kt_error, 100)


    

In [ ]: