In [ ]:
%matplotlib inline
In [ ]:
from __future__ import print_function
import numpy as np
import pylab as plt
import xraylib
import tomopy
import astra
In [ ]:
def create_phantom_tooth_orig(size_x, energy, elem_tooth, elem_implant, pixel_size, isimplant):
"""
Create phantom, share is two flowers (with 3 and 6 petals)
:param size_x: size of phantom
:param energy: energy, set in keV
:param elem_tooth: Number of the tooth chemical element
:param elem_implant: Number of the chemical element
:param pixel_size: size of one pixel, set in microns
:return: 2d array of phantom
"""
xraylib.XRayInit()
phantom = np.zeros((size_x, size_x))
sx_half = size_x / 2
sq = size_x / 14
#calculate mu
density_tooth = xraylib.ElementDensity(elem_tooth)
cross_section_tooth = xraylib.CS_Total(elem_tooth, energy)
mu_tooth = density_tooth * cross_section_tooth
density_implant = xraylib.ElementDensity(elem_implant)
cross_section_implant = xraylib.CS_Total(elem_implant, energy)
mu_implant = density_implant * cross_section_implant
#buld mesh
y, x = np.meshgrid(range(size_x), range(size_x))
xx = (x - sx_half).astype('float32')
yy = (y - sx_half).astype('float32')
r = np.sqrt(xx*xx + yy*yy)
tetta = np.arctan2(yy, xx)
#make teeth
mask_tooth = r <= sq*(1 + np.cos(2*tetta) + np.sin(2*tetta)**2)
mask_tooth += (xx*xx + yy*yy) <=(0.09*size_x)**2
mask_tooth += np.roll(mask_tooth, size_x//3, axis=0) + np.roll(mask_tooth, -size_x//3, axis=0) # make 3 teeth
phantom[mask_tooth] = mu_tooth
#make implant
mask_implant = (xx / (0.11*size_x))**2 + (yy / (0.07*size_x))**2 < 1
mask_implant *= y <= sx_half
mask_implant *= ((xx / (0.11*size_x))**2 + (((yy - 0.025*size_x) / (0.07*size_x)))**2) > 1
if(isimplant):
phantom[mask_implant] = mu_implant
phantom *= pixel_size
# print("for Ca:")
# print(mu_1*pixel_size)
# print("for Ti")
# print(mu_2*pixel_size)
return phantom
def create_phantom(size):
xraylib.XRayInit()
energy = 35.0
elem1 = xraylib.SymbolToAtomicNumber('Ca')
elem2 = xraylib.SymbolToAtomicNumber('Ti')
pixel_size = 0.04
return create_phantom_tooth_orig(size, energy, elem1, elem2, pixel_size, True)
In [ ]:
size = 512
%timeit original = create_phantom(size)
In [ ]:
original = create_phantom(size)
plt.imshow(original, cmap=plt.cm.pink)
plt.colorbar()
In [ ]:
class Spectrum(object):
"""Object to describe spectrum"""
def __init__(self, energy, intensity, label):
"""
:param energy: Energy points, keV
:param intensity: Intensity in au.
:param label: Spectrum label
"""
super(Spectrum, self).__init__()
if not len(intensity) == len(energy):
raise ValueError(
'Parameters energy and intansity mast have equal length')
self.energy = energy
self.intensity = intensity
self.label = label
In [ ]:
spectra = []
spectra.append(Spectrum(
energy= [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90],
intensity = [5.50e+16, 1.10e+16, 1.10e+16, 8.00e+15, 4.60e+15, 2.70e+15, 1.75e+15,
1.25e+15, 7.50e+14, 2.00e+15, 2.80e+15, 1.00e+15, 6.00e+14, 3.80e+14,
2.70e+14, 1.60e+14, 1.00e+14],
label = '100 kev, no filter'))
spectra.append(Spectrum(
energy= [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90],
intensity = [2.00e+15, 1.10e+16, 1.10e+16, 8.00e+15, 4.60e+15, 2.70e+15, 1.75e+15,
1.25e+15, 7.50e+14, 2.00e+15, 2.80e+15, 1.00e+15, 6.00e+14, 3.80e+14,
2.70e+14, 1.60e+14, 1.00e+14],
label = '100 kev, Al 0.5 mm'))
spectra.append(Spectrum(
energy= [10, 15, 20, 25, 30, 35, 40, 43.5, 50, 55, 60, 65, 70, 75, 80, 85, 90],
intensity = [1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+14, 2.80e+14, 7.00e+14,
1.35e+15, 7.50e+14, 2.00e+15, 2.80e+15, 1.00e+15, 6.00e+14, 3.80e+14,
2.70e+14, 1.60e+14, 1.00e+14],
label = '100 kev, Cu 0.25 mm'))
spectra.append(Spectrum(
energy= [10, 11, 13, 15, 17, 20, 23, 25, 30, 33, 35, 40, 43, 45],
intensity = [3.20e+16, 1.40e+16, 6.80e+15, 7.00e+15, 6.50e+15, 5.00e+15, 3.80e+15,
3.20e+15, 1.60e+15, 1.05e+15, 7.20e+14, 3.80e+14, 2.50e+14, 1.60e+14],
label = '50 kev, no filter'))
In [ ]:
plt.figure(figsize=(10,10))
colormap = plt.cm.gist_ncar
plt.gca().set_prop_cycle('color',[colormap(i) for i in np.linspace(0, 0.9, 4)])
for spectrum in spectra:
label = spectrum.label
plt.plot(spectrum.energy, spectrum.intensity, lw=2,label=label)
plt.hold(True)
plt.grid(True)
plt.legend()
plt.show()
In [ ]:
class Phantom(object):
"""Object to describe 2D phantom"""
def __init__(self, pixel_size=None, shape=None, label=None):
"""
:param pixel_size: Pixel size mkm
:param shape: Shape of 2D data slices
:param label: Phantom label
"""
super(Phantom, self).__init__()
self.label = label
self.shape = shape
self.pixel_size = pixel_size
self.slices = {}
xraylib.XRayInit()
def add_slice(self, phantom_slice, element_name):
if self.shape is None:
raise ValueError('Phantom shape is None, please define it')
if not phantom_slice.shape == self.shape:
raise ValueError('Phantom shape is not equal phantom_slice shape')
element_atomic_number = xraylib.SymbolToAtomicNumber(element_name)
slice_name = element_name
if not slice_name in self.slices:
self.slices[slice_name] = phantom_slice
else:
self.slices[slice_name] += phantom_slice
def get_attanuation_map(self, energy):
"""
:param energy: Energy in kev
"""
res = np.zeros(shape = self.shape)
for element_name, element_map in self.slices.iteritems():
element_atomic_number = xraylib.SymbolToAtomicNumber(element_name)
cs = xraylib.CS_Total(element_atomic_number, energy)
# print(element_atomic_number)
res += element_map*cs
return res
In [ ]:
def create_phantom_tooth(size_x, element_name='Ca'):
"""
Create phantom, share is two flowers (with 3 petals)
:param size_x: size of phantom
:param element_name: Name of the tooth chemical element
:return: 2d array of phantom
"""
xraylib.XRayInit()
phantom = np.zeros((size_x, size_x))
sx_half = size_x / 2
sq = size_x / 14
#calculate density
elem_tooth = xraylib.SymbolToAtomicNumber(element_name)
density_tooth = xraylib.ElementDensity(elem_tooth)
#buld mesh
y, x = np.meshgrid(range(size_x), range(size_x))
xx = (x - sx_half).astype('float32')
yy = (y - sx_half).astype('float32')
r = np.sqrt(xx*xx + yy*yy)
tetta = np.arctan2(yy, xx)
#make teeth
mask_tooth = r <= sq*(1 + np.cos(2*tetta) + np.sin(2*tetta)**2)
mask_tooth += (xx*xx + yy*yy) <=(0.09*size_x)**2
mask_tooth += np.roll(mask_tooth, size_x//3, axis=0) + np.roll(mask_tooth, -size_x//3, axis=0) # make 3 teeth
phantom[mask_tooth] = density_tooth
return phantom
In [ ]:
def create_phantom_implant(size_x, element_name='Ti'):
"""
Create phantom implant
:param size_x: size of phantom
:param element_name: Name of the implant chemical element
:return: 2d array of phantom
"""
xraylib.XRayInit()
phantom = np.zeros((size_x, size_x))
sx_half = size_x / 2
sq = size_x / 14
elem_implant = xraylib.SymbolToAtomicNumber(element_name)
density_implant = xraylib.ElementDensity(elem_implant)
#buld mesh
y, x = np.meshgrid(range(size_x), range(size_x))
xx = (x - sx_half).astype('float32')
yy = (y - sx_half).astype('float32')
r = np.sqrt(xx*xx + yy*yy)
tetta = np.arctan2(yy, xx)
#make implant
mask_implant = (xx / (0.11*size_x))**2 + (yy / (0.07*size_x))**2 < 1
mask_implant *= y <= sx_half
mask_implant *= ((xx / (0.11*size_x))**2 + (((yy - 0.025*size_x) / (0.07*size_x)))**2) > 1
phantom[mask_implant] = density_implant
return phantom
In [ ]:
def create_phantom_sl(size_x, element_name='Al'):
"""
Create Shepp-Logan phantom
:param size_x: size of phantom
:param element_name: Name of the implant chemical element
:return: 2d array of phantom
"""
xraylib.XRayInit()
elem = xraylib.SymbolToAtomicNumber(element_name)
density = xraylib.ElementDensity(elem)
phantom = tomopy.shepp2d(size=size_x)
phantom = np.squeeze(phantom/phantom.max())
return phantom*density
In [ ]:
size = 256
phantom = Phantom(shape=(size,size), pixel_size=10e-4)
tooth = create_phantom_tooth(size, 'Ca')
implant = create_phantom_implant(size, 'Ti')
sl = create_phantom_sl(size, 'Al')
tooth[implant>0] = 0
phantom.add_slice(tooth, 'Ca')
phantom.add_slice(implant, 'Ti')
phantom.add_slice(sl, 'Al')
In [ ]:
total = None
for k,v in phantom.slices.iteritems():
if total is None:
total = v
else:
total+=v
plt.figure()
plt.imshow(v)
plt.title(k)
plt.colorbar()
plt.show()
plt.figure()
plt.imshow(total)
plt.title('Total density')
plt.colorbar()
plt.show()
In [ ]:
for e in range(1,10,1):
p = phantom.get_attanuation_map(e)
plt.figure()
plt.imshow(p)
plt.title(e)
plt.colorbar()
plt.show()
In [ ]:
plt.figure(figsize=(10,10))
colormap = plt.cm.gist_ncar
plt.gca().set_prop_cycle('color',[colormap(i) for i in np.linspace(0, 0.9, 4)])
for e_name in phantom.slices:
element_atomic_number = xraylib.SymbolToAtomicNumber(e_name)
energy = np.arange(2,100,.1)
cs = [xraylib.CS_Total(element_atomic_number, e)*xraylib.ElementDensity(element_atomic_number) for e in energy]
cs = np.array(cs)
plt.semilogy(energy, 1.e4/cs, label = e_name)
plt.legend()
plt.grid(True)
plt.xlabel('keV')
plt.ylabel('mkm')
Интенсивность прошедшего пучка $I_1$ зависит от интенсивность падающего пучка $I_0$ так: $$I_1=I_0 \int_l {-\mu(x)dx},$$ при этом, т.к. излучение монохроматично, то $$I=N \times E,$$ где $I$ - интенсивнось излучения (поток энергии), $N$ - число упавших квантов, $E$ - энергия одного кванта.
Отсюда следует, что независимо от того, регистрируется-ли суммарная энергия или число квантов, синограмма выглядит так: $$s = I_1/I_0 = \int_l {-\mu(x)dx} $$
Если излучение полихроматичное, то: $$I_1(E)=I_0(E) \int_l {-\mu(x,E)dx},$$
Суммарное излучение на детекторе: $$I=\int_0^{E_{max}} I(E) dE $$
In [ ]:
spectr = [s for s in spectra if s.label == '100 kev, Al 0.5 mm'][0]
In [ ]:
spectr.label
In [ ]:
def create_sinogram(data, angles):
detector_size = int(data.shape[0]*np.sqrt(2))
vol_geom = astra.create_vol_geom(data.shape[0], data.shape[1])
proj_geom = astra.create_proj_geom('parallel', 1.0, detector_size, angles)
proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
W = astra.OpTomo(proj_id)
P = data
sinogram = W * P
sinogram = sinogram.reshape([len(angles), detector_size])
return sinogram
In [ ]:
for i in range(1,10):
p = phantom.get_attanuation_map(i)
s = create_sinogram(p, np.linspace(-np.pi/2,np.pi/2,180,False))
# plt.figure()
# plt.imshow(s)
# # plt.colorbar()
# plt.title(i)
# plt.savefig('{}.png'.format(i))
np.savetxt('{}.txt'.format(i),
s,fmt='%d',
header='# {} {}'.format(s.shape[0],s.shape[1]))
In [ ]: