-
Notifications
You must be signed in to change notification settings - Fork 19
/
Copy pathDoseCalculation.py
266 lines (233 loc) · 10.6 KB
/
DoseCalculation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
# -*- coding: utf-8 -*-
"""
We'd like to know a bit more about the dose we inflict on the patient.
This script is used to calculate said dose based on the x-ray spectra that we
will be able to set (see Source-Specifications).
"""
from __future__ import division # fix integer division
from optparse import OptionParser
import sys
import os
import numpy as np
from scipy import constants
# Clear command line
os.system('clear')
# Use Pythons Optionparser to define and read the options, and also
# give some help to the user
parser = OptionParser()
usage = "usage: %prog [options] arg"
parser.add_option('-v', '--kv', dest='kV',
type='float',
metavar='53',
default=90,
help='Tube peak voltage [kV] you would like to calcuate the '
'dose for. The script only accepts voltages that are '
'in the specs (and tells you if you set others). '
'Defaults to %default kV, which is the WHO setting for '
'lumbar spine.')
parser.add_option('-m', '--mas', dest='mAs',
type='float',
metavar='1.6',
default=125,
help='mAs settings. Defaults to %default mAs, which is the '
'WHO setting for lumbar spine.')
parser.add_option('-e', '--exposuretime', dest='Exposuretime',
type='float',
metavar='100',
default=1000,
help='Exposure time [ms]. Defaults to 1 second, because we '
'assume that "-m" (mAs) is used as input. If the user '
'insists, an exposure time can be set.')
parser.add_option('-d', '--distance', dest='Distance',
type='float',
metavar='100',
default=140,
help='Source-Detector distance [cm]. Defaults to %default'
'cm')
parser.add_option('-l', '--length', dest='Length',
type='float',
metavar='15',
default=43.,
help='Length of the (square) FOV [cm]. Defaults to %default '
'cm.')
parser.add_option('-t', '--thickness', dest='Thickness',
type='float',
metavar='13',
default=15.,
help='Patient or sample thickness [cm]. Used to calculate '
'attenuation. Defaults to %default cm.')
parser.add_option('-c', '--chatty', dest='chatty',
default=False, action='store_true',
help='Be chatty. Default: Tell us only the relevant stuff.')
(options, args) = parser.parse_args()
# show the help if no parameters are given
if options.kV is None:
parser.print_help()
print 'Example:'
print 'The command below calculates the dose for a peak tube voltage of', \
'60 kV.'
print
print sys.argv[0], '-v 60'
exit(1)
# Inform the user that we only have certain values to work with
Voltage = [46, 53, 60, 70, 80, 90, 100, 120]
if options.kV not in Voltage:
print 'You can only enter one of these voltages:', \
str(Voltage).strip('[]'), 'kV'
print
print 'Try again with the nearest allowed value:'
# http://stackoverflow.com/a/9706105/323100
print sys.argv[0], '-v', Voltage[min(range(len(Voltage)),
key=lambda i: abs(Voltage[i] -
options.kV))]
exit(1)
ChosenVoltage = Voltage.index(options.kV)
# Load spectra
SpectraPath = 'Spectra'
# Construct file names, then load the data with the filenames (we could do this
# in one step, but like this it's easier to debug. 'SpectrumData' is the data
# without comments, thus we read the mean energy on line 7 in a second step
SpectrumLocation = [os.path.join(SpectraPath, 'Xray-Spectrum_' +
str("%03d" % kV) + 'kV.txt')
for kV in Voltage]
SpectrumData = [(np.loadtxt(FileName)) for FileName in SpectrumLocation]
MeanEnergy = [float(open(FileName).readlines()[5].split()[3]) for FileName in
[os.path.join(SpectraPath, 'Xray-Spectrum_' + str("%03d" % kV) +
'kV.txt') for kV in Voltage]]
if options.chatty:
for v, e in zip(Voltage, MeanEnergy):
print 'Peak tube voltage', v, 'kV = mean energy', int(round(e)), 'keV'
print 'For a peak tube voltage of', options.kV, 'kV and a current of', \
int(round(options.mAs * (options.Exposuretime / 1000.))), 'mAs (exp.', \
'time', options.Exposuretime, 'ms) we get a mean energy of', \
round(MeanEnergy[ChosenVoltage], 3), 'keV.'
print
# Calculate the numbers of photons emitted from the tube.
PhotonEnergy = (MeanEnergy[ChosenVoltage] * 1000) * constants.e # Joules
# DEBUG
# to get the same value as Zhentians calculation, also change eta to 1 instead
# of 1.1, that's an error in his document :)
# options.kV = 40
# PhotonEnergy = (options.kV * 1000) * constants.e # Joules
# python DoseCalculation.py -v 46 -m 25 -e 1000 -d 120 -l 5
# DEBUG
print 'At this mean energy, a single photon has an energy of', \
'%.3e' % PhotonEnergy, 'J.'
print
# Surface entrance dose
# The K-value is based on the machine. The BAG-calculator (see below) list 0.1
K = 0.1 # mGy m^2 mAs^-1
# BSF found by Arouna2000, cited by BAG2012. Gives the same SED as the
# XLS-calculator from BAG (http://is.gd/oTpniQ)
BSF = 1.35
# calculating while converting Focusdistance from m to cm
SED = K * (options.kV / 100) ** 2 * options.mAs * (100 / options.Distance) ** 2 * BSF
print 'The surface entrance dose for an x-ray pulse with'
print '\t* U =', options.kV, 'kV'
print '\t* Q =', options.mAs, 'mAs'
print '\t* FOD =', options.Distance / 100, 'm'
print '\t* A =', int(options.Length ** 2), 'cm^2'
print '\t* K =', K, 'mGy*m^2/mAs'
print '\t* BSF =', BSF
print 'is SED = K*(U/100)^2*Q*(1/FOD)^2*BSF =', round(SED, 3), 'mGy (mJ/kg).'
print
# Correspond SED to photon count
N0 = SED / PhotonEnergy
print 'A SED of', '%.3e' % (SED / 1000), 'Gy (mJ/kg) corresponds to', \
'%.3e' % N0, 'absorbed photons per kg (with a photon', \
'energy of', '%.3e' % PhotonEnergy, 'J per photon).'
# Calculate the number of photons from the tube to the sample
# N0 = (VI/E)*eta*(A/4Pir^2)
# Calculate efficiency for a Tungsten anode according to Krestel1990, chapter
# 3.1.5
eta = 1.1e-9 * 74 * options.kV * 1000
N0 = (options.kV * 1000 * ((options.mAs / 1000) / (options.Exposuretime / 1000)) / (PhotonEnergy)) * eta * ((options.Length / 100) ** 2 / (4 * np.pi * (options.Distance / 100) ** 2))
print 'The source emits %.3e' % N0, 'photons with a mean energy of', \
'%.3e' % PhotonEnergy, 'each'
print 'We assume these photons are all the photons that reached the ' \
'patient, and thus can calculate the photon flux from this.'
Flux = N0 / (options.Exposuretime / 1000)
print 'With an exposure time of', options.Exposuretime, \
'ms the aforementioned number of photons corresponds to a photon flux ' \
'of', '%.3e' % Flux, 'photons per second (from the source to the ' \
'patient surface).'
exit()
# Attenuation in Patient
AttenuationCoefficient = 0.5 # For calculation we just simply assume 50%.
# We NEED to read the data from the NIST tables, but they're in shutdown now...
print 'Attenuation coefficient set to', AttenuationCoefficient, \
'cm^-1 (@' + str(Voltage[ChosenVoltage]), 'kV)'
# Number of absorbed photons
# N = N0(e^-uT)
N = N0 * (np.exp((-AttenuationCoefficient * (options.Thickness / 100))))
print 'Assuming an attenuation coefficient of', AttenuationCoefficient, \
'and a penetration depth of', options.Thickness, \
'cm we have (according to the Beer-Lambert law (N = N0 * e^-uT)'
print ' *', '%.3e' % N, 'photons after the xrays have passed the patient'
print ' * thus', '%.3e' % (N0 - N), 'photons were absorbed'
print ' * the intensity dropped to', round((N / N0) * 100, 2), '%'
print
print
print 'Use nist-attenuation-scraper.py to get the correct attenuation!'
# Attenuation Coefficients
# @40kV, half bone, half muscle
AttenuationCoefficient = []
AttenuationCoefficient.append(np.mean((2.685e-1, 6.655e-1)))
# @70kV (0.5*60+0.5*80), both half bone, half muscle
AttenuationCoefficient.append(np.mean((np.mean((2.048e-01, 3.148e-01)),
np.mean((1.823e-01, 2.229e-01)))))
# Skeletal muscle (http://is.gd/D88OFv)
# Energy mu/rho mu_en/rho
# (MeV) (cm2/g) (cm2/g)
# 1.00000E-02 5.356E+00 4.964E+00
# 1.50000E-02 1.693E+00 1.396E+00
# 2.00000E-02 8.205E-01 5.638E-01
# 3.00000E-02 3.783E-01 1.610E-01
# 4.00000E-02 *2.685E-01* 7.192E-02
# 5.00000E-02 2.262E-01 4.349E-02
# 6.00000E-02 *2.048E-01* 3.258E-02
# 8.00000E-02 *1.823E-01* 2.615E-02
# 1.00000E-01 1.693E-01 2.544E-02
# 1.50000E-01 1.492E-01 2.745E-02
# 2.00000E-01 1.358E-01 2.942E-02
# Cortical bone (http://is.gd/2176eQ)
# Energy mu/rho mu_en/rho
# (MeV) (cm2/g) (cm2/g)
# 1.00000E-02 2.851E+01 2.680E+01
# 1.50000E-02 9.032E+00 8.388E+00
# 2.00000E-02 4.001E+00 3.601E+00
# 3.00000E-02 1.331E+00 1.070E+00
# 4.00000E-02 *6.655E-01* 4.507E-01
# 5.00000E-02 4.242E-01 2.336E-01
# 6.00000E-02 *3.148E-01* 1.400E-01
# 8.00000E-02 *2.229E-01* 6.896E-02
# 1.00000E-01 1.855E-01 4.585E-02
# 1.50000E-01 1.480E-01 3.183E-02
# 2.00000E-01 1.309E-01 3.003E-02
r = 140 # cm, Distance from source to sample
eta = 1e-9 # *ZV
Z = 74 # Tungsten
eV = 1.602e-19 # J
QFactor = 1 # http://en.wikipedia.org/wiki/Dosimetry#Equivalent_Dose
WeightingFactor = 0.12 # http://en.wikipedia.org/wiki/Dosimetry#Effective_dose
ExposureTime = 1000e-3 # s
Current = (options.mAs / 1000) / (options.Exposuretime / 1000)
Area = (options.Length / 100) ** 2
Weight = 10 # kg
# Calculate the number of photons from the tube to the sample
# ~ N0 = (VI/E)*eta*(A/4*Pi*r^2)
N0 = (Voltage * Current) / (Voltage * eV) * eta * Z * Voltage * Area / (4 * np.pi * r ** 2)
print ' - the tube emitts %.4e' % N0, 'photons per second'
# Absorbed radiation dose per second
# Da = Eneregy / Weight # J/kg per second
Da = N * MeanEnergy * 1000 * eV / Weight
print ' -', round(Da * 1000, 4), 'mGy/s are absorbed by the sample,', \
' if we assume it is', Weight, 'kg'
# Effective dose per second
# De = Da * Wr, WR = Q * N
De = Da * QFactor * WeightingFactor
print ' -', round(De * 1000, 4), 'mSv/s is the effective dose'
# Total effective dose on the sample
D = De * ExposureTime
print ' -', round(D * 1000, 4), 'mSv is the effective dose on the', \
'sample for an exposure time of =', ExposureTime, 's)'