Simulation Models and Software Architecture - xcist/documentation GitHub Wiki

Simulation Models

A typical CatSim simulation model is given by:

where:

  • yi is the detector signal at sinogram index i,
  • k is the energy index, s is the beam sub-sampling index,
  • Aik is the number of photons arriving at the detector without attenuator for energy bin k,
  • Iiso is the intersection length between the line with index is and the object with index o,
  • μok is the linear attenuation coefficient of object o at energy k,
  • yikscatter is the scatter signal, for example computed by the Monte Carlo simulation,
  • DQE is the detector quantum efficiency (fraction of photons absorbed),
  • fCONV is a factor to convert from keV to the number of electrons, and
  • σelectronic is the standard deviation of the electronic noise.

Performing all of these requires many different models, applied at different steps in the simulation process. Here we summarize the most commonly-used simulation modes and the output files they produce.

CatSim was originally developed to produce sinograms to be used in the process of developing new reconstruction algorithms. Therefore, the focus is on simulation, rather than reconstruction. However, CatSim is also used by non-reconstruction experts, for producing reconstructed images following simulation. Therefore, some basic “prep” and reconstruction functionality is included with CatSim.

In general, these can be invoked after running CatSim, and the same configuration file can be used for prep and reconstruction as was used for the simulation itself, provided that the fields required by catprep and catsimrecon are specified in the configuration file. Details regarding configuration files are provided in the following software architecture section.

Software architecture

CatSim is written in Python (py-files) and C (C files), and defines simulations with five types of configuration files (CFG). The software architecture is shown below in Figure 1.


Figure 1: CatSim software architecture.

The intersection lengths liso described above are computed in the projector. This is CatSim’s engine and it is written in C for efficiency. In the case of analytic objects, all object intersections for a given line are sorted and the incremental path lengths are computed. For voxelized objects, path lengths are computed using the Distance Driven Projector. The Monte Carlo scatter and dose simulations also support voxelized objects.

The CT architecture, physics and acquisition are written in Python. The design is highly modular, so that – for example – different detector models can be readily applied to the spectrum present at the surface of the detector. In particular, CatSim uses the notion of "callbacks", which are user-provided scripts that can be inserted into the processing chain. These callbacks allow CatSim to simulate any acquisition geometry of interest, as well as to model any focal spots, detector geometries, etc.

The two levels (Python and C-code) are talking to each other through a clearly defined interface, which makes it straightforward to plug in new phantom types and the corresponding projector.

File types and directories

CatSim uses a number of different input and output files, which by default (but not necessarily) are located in the different sub-directories, beneath the above-listed directories.

bin/	      compiled C-code
bowtie/	      bowtie files
docs/	      documentation
materials/    material definitions
phantoms/     software phantom definitions
spectrum/     spectrum data files
mfiles/	      core CatSim mfiles – do not modify !
samples/      sample configuration scripts
src/	      C source code

The syntax of the different file types is discussed in subsequent sections.

Key configuration parameter definitions

Simulation experiments are defined by these five configuration files

  • Scanner_name.cfg: configuration file that defines the scanner itself, including the geometry and some inherent parameters of the particular scanner
  • Phantom_name.cfg: configuration file that defines the phantom to be 'scanned'
  • Protocol_name.cfg: configuration file that defines the scan protocol, including scan technique, spectrum, filtration, etc.
  • Physics_name.cfg: configuration file that defines the simulation-level parameters that determine the accuracy at which the simulation is performed
  • Recon_name.cfg: configuration file that defines the reconstruction parameters

We recommend sticking to these configuration filename conventions and parameter grouping. This will allow for easily mix-and-matching phantoms, scanner models, etc.

The specific configuration parameters include:

#* Scanner geometry scanner.detectorCallback = "Detector_ThirdgenCurved" # name of function that defines the detector shape and model scanner.sid = 540.0 # source-to-iso distance (mm) scanner.sdd = 950.0 # source-to-detector distance (mm) scanner.detectorColsPerMod = 1 # number of detector columns per module scanner.detectorRowsPerMod = 16 # number of detector rows per module scanner.detectorColOffset = 0.0 # detector column offset relative to centered position (in detector columns) scanner.detectorRowOffset = 0.0 # detector row offset relative to centered position (in detector rows) scanner.detectorColSize = 1.0 # detector column pitch or size (mm) scanner.detectorRowSize = 1.0 # detector row pitch or size (mm) scanner.detectorColCount = 900 # total number of detector columns scanner.detectorRowCount = scanner.detectorRowsPerMod # total number of detector rows scanner.detectorPrefilter = ['Al', 0.1, 'water', 2] # material file name(s) and material thickness(es) (mm)

#* X-ray tube scanner.focalspotCallback = "Source_Uniform" # name of function that defines the focal spot shape and model scanner.targetAngle = 7.0 # target angle relative to scanner XY-plane (in degrees) scanner.focalspotWidth = 1.0 scanner.focalspotLength = 1.0

#* Detector scanner.detectorMaterial = "Lumex" # detector sensor material scanner.detectorDepth = 3.0 # detector sensor depth (mm) scanner.detectionCallback = "Detection_EI" # name of function that defines the detection process (conversion from X-rays to detector signal) scanner.detectionGain = 15.0 # factor to convert energy to electrons (electrons / keV) scanner.detectorColFillFraction = 0.9 # active fraction of each detector cell in the column direction scanner.detectorRowFillFraction = 0.9 # active fraction of each detector cell in the row direction scanner.eNoise = 5000.0 # standard deviation of Gaussian electronic noise (in electrons)

#* Photon counting CT scanner.detectorMaterial = "CZT" # detector sensor material scanner.detectorDepth = 1.6 # detector sensor depth (in mm) scanner.detectionCallback = "Detection_PC" # name of function that defines the detection process (conversion from X-rays to detector signal) scanner.detectionResponseFilename = 'PC_spectral_response_CZT0.25x0.25x1.6.mat' # name of the response data file scanner.detectorBinThreshold = [20, 30, 40, 60, 80, 100, 160] # energy thresholds (keV), n bins has n+1 thresholds; the first and last are the min and max energy thresholds. scanner.detectorSumBins = 0 # 1: sum all bins (gray scale output), data dim [view row col]; 0: output multiple bins [view row col bin]

#* simulations sequence protocol.scanTypes = [1, 1, 1] # flags for airscan, offset scan, phantom scan

#* Table and gantry protocol.scanTrajectory = "Gantry_Helical" # name of the function that defines the scanning trajectory and model protocol.viewsPerRotation = 1000 # total numbers of view per rotation protocol.viewCount = 1000 # total number of views in scan protocol.startViewId = 0 # index of the first view in the scan protocol.stopViewId = protocol.startViewId+protocol.viewCount-1 # index of the last view in the scan protocol.airViewCount = 1 # number of views averaged for air scan protocol.offsetViewCount = 1 # number of views averaged for offset scan protocol.rotationTime = 1.0 # gantry rotation period (in seconds) protocol.rotationDirection = 1 # gantry rotation direction (1=CW, -1 CCW, seen from table foot-end) protocol.startAngle = 0 # relative to vertical y-axis (n degrees) protocol.tableSpeed = 0 # speed of table translation along positive z-axis (mm/sec) protocol.startZ = 0 # start z-position of table protocol.tiltAngle = 0 # gantry tilt angle towards negative z-axis (in degrees) protocol.wobbleDistance = 0.0 # focalspot wobble distance protocol.focalspotOffset = [0, 0, 0] # focalspot position offset

#* X-ray tube technique and filtration protocol.mA = 200 # tube current (in mA) protocol.spectrumCallback = "Spectrum" # name of function that reads and models the X-ray spectrum protocol.spectrumFilename = "tungsten_tar7_120_filt.dat" # name of the spectrum file protocol.spectrumScaling = 1 # scaling factor such that spectrum is in photons / mA / s / mm^2 at 1000 mm protocol.bowtie = "medium" # name of the bowtie file (or []) protocol.filterCallback = "Xray_Filter" # name of function to compute additional filtration protocol.flatFilter = ['Al',0.1] # additional filtration - materials and thicknesses (mm) protocol.dutyRatio = 1.0 # tube ON time fraction (for pulsed tubes)

#* Phantom phantom.callback = "Phantom_Voxelized" # name of function that reads and models phantom phantom.projectorCallback = "C_Projector_Voxelized" # name of function that performs projection through phantom phantom.filename = 'BIG.json' # phantom filename phantom.centerOffset = [0, 0, 0] # offset of phantom center relative to origin (mm)

#* Geometric and energy sampling physics.energyCount = 12 physics.monochromatic = -1 # monochromatic energy in kVp. if negative, means polychromatic. physics.colSampleCount = 2 physics.rowSampleCount = 2 physics.srcXSampleCount = 2 physics.srcYSampleCount = 2 physics.viewSampleCount = 2

#* Flags to determine what has to be recalculated each view physics.recalcDet = 0 physics.recalcSrc = 0 physics.recalcRayAngle = 0 physics.recalcSpec = 0 physics.recalcFilt = 0 physics.recalcFlux = 0 physics.recalcPht = 0 physics.recalcDet = 0

#* Noise on/off settings physics.enableQuantumNoise = 1 physics.enableElectronicNoise = 1

#* Internal physics models physics.rayAngleCallback = "Detector_RayAngles_2D" physics.fluxCallback = "Detection_Flux" physics.scatterCallback = "" physics.prefilterCallback = "Detection_prefilter" physics.crosstalkCallback = "" physics.lagCallback = "" physics.opticalCrosstalkCallback = "" physics.DASCallback = "Detection_DAS"

#* I/O preferences physics.outputCallback = "WriteRawView"

#* Reconstruction recon.fov = 500.0 # diameter of the reconstruction field-of-view (mm) recon.imageSize = 512 # number of columns and rows to be reconstructed (square) recon.sliceCount = 1 # number of slices to reconstruct recon.sliceThickness = 1.0 # reconstruction inter-slice interval (mm) recon.centerOffset = [0.0,0.0,0.0] # reconstruction offset relative to center of rotation (mm)

How to create new callback functions

(Future content)

Simulation Results

The output/results of CatSim are binary data in float32.

All results of CatSim and related commands are written to files for which the base name is specified in the cfg structure.

cfg.resultsName defines the basename of the result files, which can be full path and name or just the name.

The result files are:

[basename].air - air scan data, dim: [Ebin, column, row, view_air] <br>
[basename].offset - offset scan data (X-ray off), dim: [Ebin, column, row, view_offset] <br>
[basename].scan - phantom scan data, dim: [Ebin, column, row, view] <br>
[basename].prep - prep data, dim: [Ebin, column, row, view] <br>
[basename]_512x512x4.raw - reconstructed image, dim: [recon.imageSize, recon.imageSize, recon.sliceCount], here is 512x512x4 <br>
[basename].recon	reconstructed image or volume

where

  • Ebin is the energy bin, for energy integration CT, it's 1.
  • column and row are detector pixels in x and z direction.
  • view is the scan views, view_air and view_offset are 1 by default.

Basename can (optionally) include a full path. If a Basename is provided without a full path, then the Basename will be interpreted relative to the current Matlab working directory.

Output files can be viewed using the ImageJ application, described on this page.


Back to Main menu

⚠️ **GitHub.com Fallback** ⚠️