Source code for flextomo.model

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: kostenko

This module can be used to simulate spectral effects and resolution/photon-count effects.
NIST data is used (embedded in xraylib module) to simulate x-ray spectra of compounds
"""
# >>>>>>>>>>>>>>>>>>>>>>>>>>>> Imports >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

import numpy
import warnings

from flextomo import projector
from flexdata import display

# >>>>>>>>>>>>>>>>>>>>>>>>>>>> Constants >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

# Some useful physical constants:
phys_const = {'c': 299792458, 'h': 6.62606896e-34, 'h_ev': 4.13566733e-15, 'h_bar': 1.054571628e-34, 'h_bar_ev': 6.58211899e-16,
              'e': 1.602176565e-19, 'Na': 6.02214179e23, 're': 2.817940289458e-15, 'me': 9.10938215e-31, 'me_ev': 0.510998910e6}
const_unit = {'c': 'm/c', 'h': 'J*S', 'h_ev': 'e*Vs', 'h_bar': 'J*s', 'h_bar_ev': 'eV*s',
              'e': 'colomb', 'Na': '1/mol', 're': 'm', 'me': 'kg', 'me_ev': 'ev/c**2'}

# >>>>>>>>>>>>>>>>>>>>>>>>>>>> Methods >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

[docs] def ctf(shape, mode = 'gaussian', parameter = (1, 1)): """ Get a CTF (fft2(PSF)) of one of the following types: gaussian, dual_ctf, fresnel, TIE Args: shape (list): shape of a projection image mode (str): 'gaussian' (blurr), 'dual_ctf', fresnel, TIE (phase contrast) parameter (list / float): PSF/CTF parameters. For gaussian: [detector_pixel, sigma] For dual_ctf/tie: [detector_pixel, energy, src2obj, det2obj, alpha] For fresnel: [detector_pixel, energy, src2obj, det2obj] """ # Some constants... h_bar_ev = 6.58211899e-16 c= 299792458 if mode == 'gaussian': # Gaussian CTF: pixel = parameter[0] sigma = parameter[1] u = _w_space_(shape, 0, pixel) v = _w_space_(shape, 1, pixel) ctf = numpy.exp(-((u[:, None] * sigma) ** 2 + (v[None, :] * sigma) ** 2)/2) elif mode == 'dual_ctf': # Dual CTF approximation phase contrast propagator: pixelsize = parameter[0] energy = parameter[1] r1 = parameter[2] r2 = parameter[3] alpha = parameter[4] # Effective propagation distance: m = (r1 + r2) / r1 r_eff = r2 / m # Wavenumber; k = energy / (h_bar_ev * c) # Frequency square: w2 = _w2_space_(shape, pixelsize) #return -2 * numpy.cos(w2 * r_eff / (2*k)) + 2 * (alpha) * numpy.sin(w2 * r_eff / (2*k)) ctf = numpy.cos(w2 * r_eff / (2*k)) - (alpha) * numpy.sin(w2 * r_eff / (2*k)) elif mode == 'fresnel': # Fresnel propagator for phase contrast simulation: pixelsize = parameter[0] energy = parameter[1] r1 = parameter[2] r2 = parameter[3] # Effective propagation distance: m = (r1 + r2) / r1 r_eff = r2 / m # Wavenumber; k = energy / (h_bar_ev * c) # Frequency square: w2 = _w2_space_(shape, pixelsize) ctf = numpy.exp(1j * w2 * r_eff / (2*k)) elif mode == 'tie': # Transport of intensity equation approximation of phase contrast: pixelsize = parameter[0] energy = parameter[1] r1 = parameter[2] r2 = parameter[3] alpha = parameter[4] # Effective propagation distance: m = (r1 + r2) / r1 r_eff = r2 / m # Wavenumber; k = energy / (h_bar_ev * c) # Frequency square: w2 = _w2_space_(shape, pixelsize) ctf = 1 - alpha * w2 * r_eff / (2*k) return numpy.fft.fftshift(ctf)
def _w_space_(shape, dim, pixelsize): """ Generate spatial frequencies along dimension dim. """ # Image dimentions: sz = numpy.array(shape) * pixelsize # Frequency: xx = numpy.arange(0, shape[dim]) - shape[dim]//2 return 2 * numpy.pi * xx / sz[dim] def _w2_space_(shape, pixelsize): """ Generates the lambda*freq**2*R image that can be used to calculate phase propagator at distance R, photon wavelength lambda. """ # Frequency squared: u = _w_space_(shape, 0, pixelsize) v = _w_space_(shape, 1, pixelsize) return (u**2)[:, None] + (v**2)[None, :]
[docs] def apply_noise(image, mode = 'poisson', parameter = 1): """ Add noise to the data. Args: image (numpy.array): image to apply noise to mode (str): poisson or normal parameter (float): norm factor for poisson or a standard deviation """ if mode == 'poisson': return numpy.random.poisson(image * parameter) elif mode == 'normal': return numpy.random.normal(image, parameter) else: raise ValueError('Me not recognize the mode! Use normal or poisson!')
[docs] def effective_spectrum(energy = None, kv = 90, filtr = {'material':'Cu', 'density':8, 'thickness':0.1}, detector = {'material':'Si', 'density':5, 'thickness':1}): """ Generate an effective specturm of a CT scanner. """ # Energy range: if not energy: energy = numpy.linspace(10, 90, 9) # Tube: spectrum = bremsstrahlung(energy, kv) # Filter: if filtr: spectrum *= total_transmission(energy, filtr['material'], rho = filtr['density'], thickness = filtr['thickness']) # Detector: if detector: spectrum *= scintillator_efficiency(energy, detector['material'], rho = detector['density'], thickness = detector['thickness']) # Normalize: spectrum /= (energy*spectrum).sum() return energy, spectrum
[docs] def spectralize(proj, kv = 90, n_phot = 1e8, specimen = {'material':'Al', 'density': 2.7}, filtr = {'material':'Cu', 'density':8, 'thickness':0.1}, detector = {'material':'Si', 'density':5, 'thickness':1}): """ Simulate spectral data. """ # Generate spectrum: energy, spectrum = effective_spectrum(kv, filtr, detector) # Get the material refraction index: mu = linear_attenuation(energy, specimen['material'], specimen['density']) # Display: display.plot(energy, spectrum, title = 'Spectrum') display.plot(energy, mu, title = 'Linear attenuation') # Simulate intensity images: counts = numpy.zeros_like(proj) for ii in range(len(energy)): # Monochromatic component: monochrome = spectrum[ii] * numpy.exp(-proj * mu[ii]) monochrome = apply_noise(monochrome, 'poisson', n_phot) / n_phot # Detector response is assumed to be proportional to E counts += energy[ii] * monochrome return counts
[docs] def forward_spectral(vol, proj, geometry, materials, energy, spectrum, n_phot = 1e8): """ Simulate spectral data using labeled volume. """ max_label = int(vol.max()) if max_label != len(materials): raise ValueError('Number of materials is not the same as the number of labels in the volume!') # Normalize spectrum: spectrum /= (spectrum * energy).sum() # Simulate intensity images: lab_proj = [] for jj in range(max_label): # Forward project: proj_j = numpy.zeros_like(proj) vol_j = numpy.float32(vol == (jj+1)) projector.forwardproject(proj_j, vol_j, geometry) lab_proj.append(proj_j) for ii in range(len(energy)): # Monochromatic components: monochrome = numpy.ones_like(proj) for jj in range(max_label): mu = linear_attenuation(energy, materials[jj]['material'], materials[jj]['density']) monochrome *= numpy.exp(-lab_proj[jj] * mu[ii]) monochrome *= spectrum[ii] monochrome = apply_noise(monochrome, 'poisson', n_phot) / n_phot # Detector response is assumed to be proportional to E proj += energy[ii] * monochrome
[docs] def material_refraction(energy, compound, rho): """ Calculate complex refrative index of the material taking into account it's density. Args: compound (str): compound chemical formula rho (float): density in g / cm3 energy (numpy.array): energy in KeV Returns: float: refraction index in [1/mm] """ import xraylib cmp = xraylib.CompoundParser(compound) # Compute ration of Z and A: z = (numpy.array(cmp['Elements'])) a = [xraylib.AtomicWeight(x) for x in cmp['Elements']] za = ((z / a) * numpy.array(cmp['massFractions'])).sum() # Electron density of the material: Na = phys_const['Na'] rho_e = rho * za * Na # Attenuation: mu = mass_attenuation(energy, compound) # Phase: wavelength = 2 * numpy.pi * \ (phys_const['h_bar_ev'] * phys_const['c']) / energy * 10 # TODO: check this against phantoms.m: phi = rho_e * phys_const['re'] * wavelength # Refraction index (per mm) return rho * (mu / 2 - 1j * phi) / 10
[docs] def mass_attenuation(energy, compound): ''' Total X-ray absorption for a given compound in cm2g. Energy is given in KeV ''' import xraylib # xraylib might complain about types: energy = numpy.double(energy) if numpy.size(energy) == 1: return xraylib.CS_Total_CP(compound, energy) else: return numpy.array([xraylib.CS_Total_CP(compound, e) for e in energy])
[docs] def linear_attenuation(energy, compound, rho): ''' Total X-ray absorption for a given compound in 1/mm. Energy is given in KeV ''' # xraylib might complain about types: energy = numpy.double(energy) # unit: [1/mm] return rho * mass_attenuation(energy, compound) / 10
[docs] def compton(energy, compound): ''' Compton scaterring crossection for a given compound in cm2g. Energy is given in KeV ''' import xraylib # xraylib might complain about types: energy = numpy.double(energy) if numpy.size(energy) == 1: return xraylib.CS_Compt_CP(compound, energy) else: return numpy.array([xraylib.CS_Compt_CP(compound, e) for e in energy])
[docs] def rayleigh(energy, compound): ''' Compton scaterring crossection for a given compound in cm2g. Energy is given in KeV ''' import xraylib # xraylib might complain about types: energy = numpy.double(energy) if numpy.size(energy) == 1: return xraylib.CS_Rayl_CP(compound, energy) else: return numpy.array([xraylib.CS_Rayl_CP(compound, e) for e in energy])
[docs] def photoelectric(energy, compound): ''' Photoelectric effect for a given compound in cm2g. Energy is given in KeV ''' import xraylib # xraylib might complain about types: energy = numpy.double(energy) if numpy.size(energy) == 1: return xraylib.CS_Photo_CP(compound, energy) else: return numpy.array([xraylib.CS_Photo_CP(compound, e) for e in energy])
[docs] def scintillator_efficiency(energy, compound='BaFBr', rho=5, thickness=1): ''' Generate QDE of a detector (scintillator). Units: KeV, g/cm3, mm. ''' # Attenuation by the photoelectric effect: spectrum = 1 - numpy.exp(- thickness * rho * photoelectric(energy, compound) / 10) # Normalize: return spectrum / spectrum.max()
[docs] def total_transmission(energy, compound, rho, thickness): ''' Compute fraction of x-rays transmitted through the filter. Units: KeV, g/cm3, mm. ''' # Attenuation by the photoelectric effect: return numpy.exp(-linear_attenuation(energy, compound, rho) * thickness)
[docs] def bremsstrahlung(energy, energy_max): ''' Simple bremstrahlung model (Kramer formula). Emax ''' if energy_max < 10: warnings.warn('Maximum energy in the Kramer formula is too low. Setting to 100kV') energy_max = 100 # Kramer: spectrum = energy * (energy_max - energy) spectrum[spectrum < 0] = 0 # Normalize: return spectrum / spectrum.max()
[docs] def gaussian_spectrum(energy, energy_mean, energy_sigma): ''' Generates gaussian-like spectrum with given mean and STD. ''' return numpy.exp(-(energy - energy_mean)**2 / (2 * energy_sigma**2))
[docs] def nist_names(): ''' Get a list of registered compound names understood by nist ''' import xraylib return xraylib.GetCompoundDataNISTList()
[docs] def find_nist_name(compound_name): ''' Get physical properties of one of the compounds registered in nist database ''' import xraylib return xraylib.GetCompoundDataNISTByName(compound_name)
[docs] def parse_compound(compund): ''' Parse chemical formula ''' import xraylib return xraylib.CompoundParser(compund)