[18JUL23 text]
Materials are fundamental to several aspects of X-ray/CT imaging systems, including source and detector filtration, the detector's X-ray conversion materials, and the patient (phantom).

In XCIST, materials are specified in a simple text format. We have provided the CatSim material calculator to produce material files in a consistent and accurate manner. Using this tool, you can quickly specify new materials based on chemical formulas or mass fractions. Detailed instructions are provided within the tool. If you need a new material, please download & use the tool, and, if not a proprietary material, please contribute your new material specification back to XCIST!

At the time of this writing (July 2023), the material folder contained 92 elements, 29 compounds and mixtures, and 40 biological tissues.

[Original text]
The ASCII text-files for the material composition and constants are located in a directory catsim/materials/ and have no extension (files with extension .mac and .rho are obsolete material definition files, and have been moved to the obsolete/materials directory). The first line specifies the number of elemental components. The second line contains the mass density (g/cm³), and the next lines contain the atomic number and the mass fraction of the respective components. The mass density is absolute. The mass fractions are relative (they don’t need to sum to 1 – they are normalized in the code).

Disclaimer: Material files included with CatSim may or may not be precisely correct. In general, you can probably assume that elemental material files are correct, but material files for mixtures may include approximations.

The mass attenuation coefficients of the materials are determined from files that are stored in subdirectories under catsim/materials. These files separate the energy-dependent Compton, photoelectric, and Rayleigh absorption processes for each element. These files were obtained from Geant4, and are used under the Geant4 Software License, found at https://geant4.web.cern.ch/download/license.

Any of the phantoms can be converted to a voxelized phantom using the function catvoxel.m. This will produce output that can be used by catdose.m and catscatter.m as well as by catsim.m.

The routine that reads and parses the phantom is specified using cfg.callback_phantom = 'PhantomCallback';
where 'PhantomCallback' is one of

'Phantom_Analytic'
'Phantom_NCAT'
'Phantom_Voxelized'

The phantom file is given by

cfg.phantom_filename = 'FileName';

where FileName is a file with extension .pp (FORBILD syntax), .ppm (CatSim syntax), or .nrb (NCAT syntax). These can be located in the phantoms directory, but a full path can be included optionally, or the file can be specified relative to the current working directory.

cfg.phantom_scale = ScaleFactor;


is an optional scale factor to be applied to the phantom. Generally CatSim uses mm for all dimensions.

There are several different phantom types/formats, including

• voxelized,
• analytic,
• polygonal, and
• NURBS (the format for files produced by XCAT software).

Each type requires its own projector because, as a ray is traced from the source to the detector, we need to determine the pathlength though each material in the phantom, as well as the volume fraction of that material. Determination of the pathlength and material properties depends on the phantom format, hence the need for projectors that can parse each format.

Voxelized files are specified as one or more 2D images. For details about the format, see

• Phantom files: Voxelized format

Because voxelized phantoms are specified as images, they can be visualized using an image viewer. For details about visualization using ImageJ, see Viewing image files using ImageJ.

The latter three phantom types are specified in text files, each with their own unique format. For details about those formats, see

• Phantom files: Analytic format
• Phantom files: Polygonal format
• Phantom files: NURBS format

These phantom types are difficult to visualize, because it would be necessary to parse the file and render an image of what it describes. At this time, no such utilities are provided with XCIST.

As explained in section 6: Software architecture, CatSim is in principle able to use different phantom types, provided a C library is provided to perform the corresponding projection operations. The projector is specified by the field

cfg.callback_projector = 'ProjectorCallbackFunctionName';

Where ProjectorCallbackFunctionName is one of

'C_Projector_Analytic'
'C_Projector_NCAT'
'C_Projector_Polygon'
'C_Projector_Voxelized'

These all have been reasonable well tested, with the exception of the polygon projector. Prior to Release 5.7, the voxelized projector was limited to a monochromatic spectrum and a single material (water). With Release 5.7, we extended the voxelized projector to perform polychromatic simulations of multi-material voxelized phantoms.

The core catsim.m routine has only a few calls to the C libraries. The various components in the projection process are parsed and passed to C in the following sequence, via the routines specified by these variables:

cfg.callback_getphantom
cfg.callback_setMaterial
cfg.callback_setSource
cfg.callback_setDetector
cfg.callback_projector

The type of phantom that is used determines the projector callback required, which in turn determines the callback that will be used for each of the others. This relationship is summarized in Table 1: Phantom types and associated callback functions.

Some projectors also use the following fields:

The maximum number of threads is specified using
cfg.n_threads=1;

The detector bounding box computation is specified using
cfg.detector_bounding_box_strategy=2;

where 1 means that only v limits are computed for the detector bounding box, and 2 means that both u and v limits are used.

Table 1: Phantom types and associated callback functions

Obviously, the first thing the voxelized projector needs is a voxelized phantom. Please read the CatVoxel section of the Wiki.

The features in an analytic phantom are continuous (not approximated as discrete linear elements), therefore these produce “ideal” simulation results – there is no relationship between the intrinsic phantom characteristics and the pixel size of the detector, or the voxel size in the reconstruction matrix. In contrast, a voxelized phantom inherently contains linear segments that might represent continuous curves. Therefore, one needs to be cautious about the precision of the voxelized phantom uses, depending on the application.

In general, we recommend that the voxelized phantom contain voxels where

• The phantom voxel size in X and Z is smaller than the detector size in X and Z, divided by the geometric magnification.
• The phantom voxel size in X, Y, and Z is smaller than the reconstruction matrix size in the same directions.

In some cases, a voxelized phantom with somewhat large voxels (e.g. > 1 mm) may be perfectly adequate for the application. But typically, the voxelized phantom should have sub-millimeter voxels.

The effect of the choice of phantom voxel size relative to reconstruction voxel size is illustrated in section Voxelized projector validation.

Because voxelized phantoms are image volumes for each material, the total memory required for a phantom of many materials could be large, especially if the phantom represents a large object (in XY and Z), is made with high resolution (many small voxels), and includes many materials.

Required memory can be calculated by

M = NumVoxelsX * NumVoxelsY * NumVoxelsZ * NumMaterials * 4 bytes per voxel

As extreme examples, a (200 mm)³ phantom with 2 materials, represented as a (100 pixel)³ volume (i.e. 2 mm voxels), will require only 8 MB of memory.

However, a (500 mm)³ phantom with 20 materials, represented as a (1k pixel)³ volume (i.e. 0.5 mm voxels), will require 80 GB of memory.

If the required memory size can fit in available RAM, then the entire phantom can be loaded at once, and the simulation will be most efficient. Alternatively, one material can be loaded into RAM at a time. However, this needs to be repeated for each view, therefore efficiency and simulation performance suffer dramatically. The behavior of the voxelized projector in this regard is specified by

cfg.enough_memory_for_voxelized_phantom = Flag;

If the flag is set to 1, CatSim will attempt to load the entire phantom, and a fatal error will occur if there is insufficient memory. In that case, the flag must be set to 0 in order to perform the simulation.

To validate the voxelized projector when using multiple materials and a polychromatic spectrum, we ran two experiments, using multiple variations in each. One experiment used an analytic phantom and the analytic and voxelized projectors; the second experiment used an NCAT phantom and the NCAT and voxelized projectors. In each experiment, we compared the results of the two projectors after reconstructing images using identical reconstruction parameters.

For the configuration files used for these experiments, see

/base/samples/test_uniformity/readme.txt
/base/samples/test_NCAT/readme.txt

All experiments were performed under the following conditions:

1. Multi-material phantoms were used.
2. A 120 kVp spectrum was applied, using (12) 10 keV bins.
3. “LightSpeed VCT” system and detector geometry were used. For some of Experiment 2, the detector pixels were binned 4X and the number of views was reduced 4X.
4. No source or detector oversampling was used for the analytic and NCAT projectors.
5. Oversampling of 11 was used when voxelizing the phantoms.
   

CatSim was originally developed to simulate the following analytic objects: cylinders, ellipsoids, boxes, and clipping planes. Subsequently, the capability to simulate other analytic objects has been added including hyperboloids of one plane, hyperboloids of two planes, cones, and tori (see figure to the right). However, the FORBILD phantom model does not support all these objects. For this and other reasons, we developed a new phantom format that’s based on Matlab syntax. We use the extension .ppm for the new (Matlab-based) format.

This syntax supports phantom files that contain one of two classes of shapes: analytic and polygonal.

For analytic phantoms, use cfg.callback_projector = 'C_Projector_Analytic';

For polygonal phantoms, use cfg.callback_projector = 'C_Projector_Polygon';

Files in the .ppm file syntax are read in by executing them as Matlab scripts. At the end of execution, it is expected that each object in the phantom will be stored in a structure called object, and a list of the materials will be stored in a cell array called materialList. The object structure contains various fields as described below.

These fields are relevant to analytic and polygonal phantoms:

These fields are relevant only to polygonal phantoms:

These fields are relevant only to analytic phantoms:

There are many ways to structure the phantom description file since it is executed within Matlab. The most straightforward (and verbose) way to specify a given object is to set each parameter of each object on a separate line. For example:

object{1}.center = [0.000000 0.000000 0.000000];
object{1}.half_axes = [108.000000 108.000000 17.500000];
object{1}.euler_angs = [0.000000 0.000000 0.000000];
object{1}.density = 1.000000;
object{1}.type = 'Cylinder';
object{1}.material = 'water';
object{1}.axial_lims = [0.000000 0.000000];
object{1}.shape = 0.000000;
object{1}.clip = [];
object{1}.color = [1.000000 0.000000 0.000000];
object{1}.group = 0;
object{1}.transparancy = 0.000000;   

obj=[];object=[];
update = '[obj,object]=AddObject(obj,object,materialList);'

If this header is included, you can specify an object in a variable called obj and then evaluate the update step to add it to the phantom (which is stored in object).

Furthermore, loops can be used to add many objects (e.g. resolution bars) with a small amount of code.

Body phantoms are described in the NCAT NURB syntax, while lesion phantoms are described as mesh, triangular polygons. Several lesion phantoms can be described in a single file.

The polygonal mesh syntax is described as follows:

• The file name starts with "tmod_" to indicate a mesh phantom
• For each lesion model described in the file:
  • The first line describes the phantom object
  • The second line denotes the material (see Table 2 below)
  • The third line is an integer describing the total number of triangles that make up the phantom
  • Each line that follows consists of three sets of x,y,z coordinates that describe the three vertices of a triangle in three-dimensional space.

Example: the file "tmod_lung_lesion_8mm.nrb" describes six different lesion phantoms made up of muscle material: 8 mm-, 7 mm-, 6 mm-, 11 mm-, another 11 mm-, and another 6 mm-long lesion. The first six lines for the 8 mm and 7 mm phantoms are excerpted from "tmod_lung_lesion_8mm.nrb" below: the 8 mm is made up of 80,636 triangles, the 7 mm is made of 49,938 triangles:

tmod_lung_lesion_8mm
1
80636
-11.842316,12.247543,122.757446 -11.817810,12.261429,122.778320 -11.836731,12.230606,122.772827 
-11.836731,12.230606,122.772827 -11.817810,12.261429,122.778320 -11.845520,12.235504,122.801025 
-11.817810,12.261429,122.778320 -11.842712,12.252777,122.817139 -11.845520,12.235504,122.801025 
...

tmod_lung_lesion_7mm
1
49938
-69.139404,8.895981,126.426636 -69.155350,8.924042,126.414673 -69.183807,8.918701,126.470337 
-69.183807,8.918701,126.470337 -69.205612,8.875168,126.515503 -69.139404,8.895981,126.426636 
-69.155350,8.924042,126.414673 -69.194672,8.933533,126.446289 -69.183807,8.918701,126.470337 
...

For these phantoms, use
cfg.callback_projector = 'C_Projector_NCAT';

Two programs, which are not part of CatSim, are used to produce phantoms in this format. These programs, developed by Paul Segars, et. al., are called NCAT and XCAT. Therefore, the resulting phantoms, and the functions associated with using them, often contain NCAT or XCAT in their names. Both NCAT and XCAT produce phantoms in the NURB format.

Note to GE’s academic collaborators: To obtain the software required to produce NCAT phantoms, please contact Paul Segars at [email protected].

NCAT and XCAT phantoms are described in files with the .nrb extension. This format is used to specify non-uniform rational b-spline (NURBS) objects for use with the ncat projector. There are m × n control points that specify a mesh that controls the shape of the NURBS object. For more information on NURBS, Wikipedia is a good resource. The files have a series of objects, with each object specified as follows:

**object_name
material_id
m :M
n :N

U Knot Vector
u(1)
...
u(n+4)

V Knot Vector
v(1)
...
v(m+4)

Control Points
c(1,1)
...
c(1,n)

c(2,1)
...
c(2,n)
...
...
c(m,1)
...
c(m,n)

The material_id is a number from 0 to 34 that specifies the material. The ncat materials are as follows:

We have implemented a "list" configuration for simulating a Hybrid projector.

An example of a hybrid projector simulation script file Sim_Sample_Hybrid.py that incorporates voxelized and analytic projectors can be found at:

(https://github.com/xcist/main/blob/master/gecatsim/examples/Sim_Sample_Hybrid.py)

An example of the corresponding configuration file Phantom_Sample_Hybrid.cfg can be found at:

(https://github.com/xcist/main/blob/master/gecatsim/examples/cfg/Phantom_Sample_Hybrid.cfg)

CatSim now provides a very basic means of producing dynamic simulations. This is done by specifying a series of phantom files, each of which provides a snapshot of the phantom at a particular time. The filenames of these phantoms must all be identical except for the specification of the phase number. For example, if simulation of a phantom with 200 unique phases is desired, 200 different files would be needed. These could, for example, be named my_dynamic_1.ppm, my_dynamic_2.ppm, my_dynamic_3.ppm, …, my_dynamic_200.ppm. In the .cfg file, the phantom filename would be specified as my_dynamic_&.ppm in this case. Other changes to the config file might be as follows:

cfg.dynamic_phantom = 1;
Nphases = 200;
bpm = 60;
cfg.phase_length = (60/bpm)*(ones(1,Nphases)/Nphases);
start_time=0;
% optional dynamic phantom parameters
% (do not specify one but not the other)
cfg.phase_start_time = start_time;
cfg.current_phase = 1;

Note that there are two required variables. The first is the cfg.dynamic_phantom flag, which must be set to 1. The other is cfg.phase_length, which is a vector that specifies the amount of time that the phantom is in each state. There are also two optional parameters (cfg.phase_start_time and cfg.current_phase) that allow one to specify a phase other than the first phase as the initial phase.

There is a second means for specifying a dynamic phantom that requires only one filename. To use this second dynamic phantom mode, you must set the following:

cfg.dynamic_phantom = 2;

This second dynamic mode is only compatible with the .ppm format. Each time the phantom is reloaded (i.e., you are in a new phase), the .ppm file is re-executed. The mechanism for making the phantom different is that the cfg parameter time is available when the .ppm file is executed. Thus, to make a sphere that moves in z sinusoidally at a rate of one period per second and an amplitude of 10 mm, you could do this:

obj.center = [0 0 10*sin(cfg.time*2*pi)];
obj.type = ‘Ellipsoid’;
obj.half_axes = [50 50 50];
eval(update);

This assumes you are using the .ppm header specified in example.ppm

If you have phantom files that are large, you will need to add it to the GitHub LFS (Large File System). First, install LFS on your computer by following these instructions: https://docs.github.com/en/repositories/working-with-files/managing-large-files/installing-git-large-file-storage?platform=windows After you've verified that LFS is installed correctly, run the following commands:

git lfs track "*.nrb"

This is to tell git which files will be stored in the LFS. In the above example, all files with the ".nrb" extension will be stored in the LFS. If your phantom files have a different extension, use that instead of ".nrb".

Next, run this command to add a phantom file. (This is similar to adding a regular file to git.)

git add <my_phantom_file_name.nrb>

You can also add multiple files at the same time as in the case of regular files:

git add *.nrb

However, because phantom files are typically quite large (50-70+ MB), it is recommended that you add one phantom file at first to make sure that it is correctly uploaded to the GitHub LFS without any issues. Once that file has been successfully uploaded, you can upload the rest of them together.

Commit your files and push, just as with regular files:

git commit -a
git push

Verify that the files you uploaded are in the LFS. There are two ways to do this:

• Go to the URL of the repository and folder via your regular web browser. Look for the file size and also a line indicating 'Stored with Git LFS.'
• From the command line, you can run this command to show all your files that in the LFS:

git lfs ls-files

The focal spot is modeled by a number of sample points, each defined by a location (x,y,z) and its relative weight w. These locations and weights are defined by a source callback function, specified by

cfg.callback_source = 'SourceCallbackFunctionName';

CatSim includes callback functions Source_Uniform.m to define a uniform (rectangular) focal spot, Source_Gaussian.m to define a Gaussian focal spot, and Source_Performix.m to define a typical bi-modal Performix focal spot. Input arguments these callback functions include:

cfg.target_angle = SourceTargetAngle_Degrees;

the angle of the target (degrees) relative to the xy-plane: a positive angle means that the plane of the focal spot is tilted away from the positive y-axis in the negative z-direction (similar to gantry tilt angle).

cfg.fs_width = FocalSpotWidth_Millimeters;

the width of the focal spot (x-direction) (mm).

cfg.fs_length = FocalSpotOpticalLength_Millimeters;

the (optical) length of the focal spot (z-direction) (mm), as viewed from the detector.

Instead of defining the (optical) length, one can also specify

cfg.fs_thermal_length = FocalSpotThermalLength_Millimeters;

the thermal length of the focal spot (yz-direction) (mm). If both are specified, the thermal length will be ignored and computed from the (optical) length.

In addition to the source locations and weights, the source callback functions also associate a planar convex hull with the focal spot, defined by a number of corner points. It is the responsibility of the user to ensure that all sample points are inside the hull, as seen from the detector. The source callback functions also associate three unit vectors with the source: front represents the central direction of the x-ray beam, lateral represents the lateral or fan direction, and long represents the longitudinal or cone direction. These are used to define a direction-dependent flux profile, source collimation, bowtie, ...

The source and focal spot are defined as a structure src containing the following fields:

src.n_samples	number of source samples
src.samples	[3 x n_samples] array of (x,y,z)
src.weights	[n_samples] array of weights
src.front	3D vector defining the central direction of the x-ray beam
src.lateral	3D vector defining the lateral (fan) direction of the x-ray beam
src.long	3D vector defining the longitudinal (cone) direction of the x-ray beam
src.n_corners	number of corners of a convex and planar focal spot region
src.corners	[3 x n_corners] array of (x,y,z) coordinates of the corners

Note: In order to take full advantage of the focal spot and detector models, and correctly model spatial resolution, appropriate oversampling must be utilized. See the Sampling section of the Wiki.

*** Also refer to validation papers ***

cfg.callback_spectrum = 'Spectrum';

specifies the function that reads or computes the spectrum. cfg.spectrum_filename specifies the file that contains the spectral information of the source. More details are given below. cfg.energy_weighting is a flag to switch between photon counting (0) and energy integration (1).

The energy dependence of the X-ray source is defined by a spectrum file. The detector dependence is defined by the scintillator material, scintillator thickness and the energy weighting. The spectrum files are ASCII text files. The first row contains an integer defining the number of energy bins in the spectrum file. All other rows contain a pair of numbers: the energy (in keV) and the number of photons in that energy bin. The first row represents the lowest energy in the spectrum, and the last row represents the highest energy in the spectrum. If needed, the spectrum will be re-binned to the desired number of energy bins. The X-ray takeoff angle may be provided in the end of the file, otherwise it always assumes that the target anode angle is its takeoff angle.

In 2011, the existing CatSim x-ray spectrum files were carefully evaluated related to the contrast and noise that were produced by equivalent simulations and measurements in uniform circular phantoms. Results, particularly related to noise, we not as closely matched as we had hoped. New spectra were then developed and evaluated using three different spectrum generation tools, at the four common kVs (80, 100, 120, 140). The spectra were also evaluated by performing simulations and measurements of the signal detected after transmission through various thicknesses of PMMA and also aluminum, using the same four spectra. Validation results are recorded in separate documentation. We evaluated the four common kVps with a 20 cm water phantom (the “QA” phantom), and 35 cm and 48 cm LDPE phantoms. We found disagreement of <5% in noise over the entire range of conditions, except that the disagreement was even greater for 80 kVp in the larger phantoms. This was a substantial improvement over the existing spectra.

The results led us to a new strategy for CatSim spectra:

1. For most accurate results, use the new spectrum files in the “spectrum” directory. These include the attenuation of the inherent filtration of the VCT/HD (“Hercules”) X-ray tube. They are named in the format xspect_tar7_kVp_filt.dat, which indicates that they were produced using XSPECT (v3.5), a 7° tungsten anode angle, kVp available from 60 to 160 kVp in 10 kVp steps, and the inherent filtration was applied. These spectra were only evaluated in the central few planes; expect some modest errors at the extreme ends of the VCT/HD cone beam.



When using these files, if you want to use the spectrum of the VCT/HD, be sure that the cfg.flat_filters parameter is not specified in your configuration file. (This will default to no added filtration.) But of course, if you want to simulate a system that has added filtration, cfg.flat_filters can be used to do that.


The older spectrum files (sp_kVp.dat and VCTkVpN.dat, as well as a few miscellaneous others) were moved to a “obsolete/spectrum” directory. If you would like to use these files, (e.g. for compatibility with older experiments) you will need to revise your configuration file to point to the “obsolete/spectrum” directory. Please do not move them into the primary “spectrum” directory in your local copy of CatSim, because then if you inadvertently “commit” these back into the CatSim code repository, they will be in the wrong place.


2. If you want to perform simulations with other conditions (e.g. other anode angles, inherent filtration, etc.), there are a few options.
   

   Bowtie and flat filters

   Bowtie filters are specified by the thickness (in cm) of Al, graphite, Cu, and Ti, as a function of alpha (the fan-angle of a detector column relative to the center channels).
   

   cfg.bowtie_filename = 'BowtieFilename';
   

   specifies the file that defines the bowtie.
   

   Bowtie definition files are text files located in the directory bowtie. Each file contains one row for each alpha, and each row contains 5 numbers: alpha, Al thickness (cm), graphite thickness (cm), Cu thickness (cm), and Ti thickness (mm). VCTlarge.txt and VCTsmall.txt represent approximations of the GE LightSpeed bowtie filters. If no bowtie is needed, specify none.txt.
   

   cfg.callback_bowtie = 'Bowtie';
   

   specifies the function that gets the bowtie definition from BowtieFilename and converts into the format desired by CatSim. GetBowtie does not support variations in z such as the Heel effect.
   

   cfg.enable_dynamic_bowtie is a flag that enables dynamic bowties.
   

   For most accurate results, the “VCT” bowties should be used. The older bowtie files (bodybow.txt, headbow.txt, and none.txt) have been moved to the obsolete/bowtie directory. If you would like to use these files, (e.g. for compatibility with older experiments) you will need to revise your configuration file to point to the obsolete/bowtie directory. Please do not move them into another bowtie directory in your local copy of CatSim, because then if you inadvertently svn commit these back into the CatSim code repository, they will be in the wrong place.
   

   Note: The bowtie file does not specify a “solid model” of the bowtie – it specifies the path length through the different materials of a bowtie. The only bowtie files implemented are for the LightSpeed VCT. For the 750HD, the bowtie is the same physical design, but it’s located at a slightly different position with respect to the source focal spot. Therefore, the bowtie files for the 750HD would have to be different, but they have not been updated to reflect this change.
   

   Flat filters can also be placed between the source and the subject. cfg.flat_filters can specify a series of materials and thicknesses to define flat pre-patient filters. The material of the filter has to correspond to material files in materials/. See the Materials section of the Wiki. For example:
   

   cfg.callback_filtration = 'Filter_Flat';
   cfg.flat_filters = {'al' 1.5 'cu' 0.075};
   

   
   specifies a filter stack of 1.5 mm aluminum and 75 microns of copper.
   
   Note: CatSim applies a cosine factor to compute the path length through these filters across alpha (in-plane angle). Because the filters are flat, there is no dependence on the filter position with respect to the source focal spot. CatSim does not apply a cosine factor across gamma (angle in the cone beam (z) direction), unless the spectrum is also defined as a function of cone angle, in which case a cosine factor in the gamma direction is also applied.

   CatVoxel

   CatVoxel will be a future feature.

   Detector and detector modules

   The detector is defined by a callback function. In principle, there can be multiple callback functions available, each with its own features and configuration variables. In practice, there is currently only one released callback function, which is specified by the following line in CatSim’s configuration file.
   

   cfg.callback_detector = 'Detector_ThirdgenCurved';
   

   This computes the entire detector structure, including the module definitions, module coordinates, cell coordinates, sample coordinates and weights. This routine uses the parameters discussed below.

   Third generation/Curved CT detector

   The detector is built as a set of planar detector modules, which can but do not have to represent actual (manufactured) detector modules. These modules are arranged side by side in an arc such that columns in the detector modules are parallel, and rows in the detector modules form a piecewise linear arc. The radius of the arc (in mm), from the nominal center of the source to the center of each detector module, is given by
   

   cfg.sdd = Source-to-Detector Distance;
   

   Only one row of modules is permitted (i.e. modules can be replicated laterally – in the row direction – but not longitudinally). The number of modules in the detector, i.e. replicates in the row direction, is specified by
   

   cfg.n_modules = number of detector modules;
   

   There will be typically only one module type, but many modules. However, it is possible to model a planar detector by specifying a single module. One disadvantage to this is that some of the projectors are multi-threaded. This is accomplished by distributing the projection task on a per-module basis. Therefore, the performance improvement of multi-threading would not be realized if a single detector module were specified.
   

   The numbers of columns and rows per module are specified by
   

   cfg.cols_per_mod = number of columns per module;
cfg.rows_per_mod = number of rows per module;
   

   The total number of rows in the detector is given by
   cfg.row_count = cfg.rows_per_mod;
   

   The total number of columns in the detector (sum for all modules) is given by
   cfg.col_count = number of detector columns;
   

   If this is specified as a constant that is less than cfg.cols_per_mod*cfg.n_modules, then the size of the “edge” modules is reduced accordingly.
   

   The total number of cells in the detector (sum for all modules) is given by
   cfg.total_n_cells = cfg.col_count*cfg.row_count;
   

   The detector column size and row size, or the u and v dimensions of a single detector cell (in mm) are specified by
   cfg.col_size = detector pixel u dimension (mm);
cfg.row_size = detector pixel v dimension (mm);
   

   The detector fill fraction is specified by
   cfg.col_fillfraction = detector pixel fill fraction in u direction;
cfg.row_fillfraction = detector pixel fill fraction in v direction;
   

   This models the fact that there is a small dead zone between detector columns and rows respectively. A fill fraction of 1.0 means that the detector has no dead space. A fill fraction of 0.0 means that the detector is completely inactive.
   

   Note: In order to take full advantage of the focal spot and detector models, and correctly model spatial resolution, appropriate oversampling must be utilized. See section 23: Sampling.
   

   The number of samples to use per detector cell in the column and row direction is specified by
   cfg.col_oversample = detector pixel oversampling in u direction;
cfg.row_oversample = detector pixel oversampling in v direction;
   

   The sample locations are distributed uniformly inside the active area of the detector. Some additional samples are placed at the edge of the neighboring detector cells to model x-ray crosstalk. X-ray crosstalk (not including optical crosstalk) is specified by
   cfg.col_crosstalk = detector pixel crosstalk fraction in u direction;
cfg.row_crosstalk = detector pixel crosstalk fraction in v direction;
   

   which is the fraction of the signal that is detected as cross-talk in neighboring detector columns/rows. For example 0.05 means that each column/row receives 5% of the x-rays from both neighboring columns/rows.
   

   If a CatSim user defines his own detector_callback routine, some additional parameters may be required.
   

   In the detector callback function, the detector is defined as a structure det containing the following fields: (see Detector definition figure below)

   
   

   det.n_modules	number of detector modules
   det.modtypes	vector of size [det.n_modules] specifying the module types; this is one-based in Matlab and will automatically be converted to zero-based when passing to C. So if all modules are of the same type, this should be an [det.n_modules] array of all ones.
   det.modcoords	[3 * det.n_modules] array of module centers (x,y,z) in mm.
   det.uvecs	[3 * det.n_modules] array of unit vectors defining a u-direction (x,y,z) for each module. The u-vector is typically perpendicular to the rotation axis.
   det.vvecs	[3 * det.n_modules] array of unit vectors defining a v-direction (x,y,z) for each module. The v-vector is typically parallel to the rotation axis.
   det.total_n_cells	total number of detector cells, i.e., the sum of the number of cells in all modules.
   det.n_moddefs	number of different module definitions; often this will be one, assuming all modules are of the same type.
   det.n_cells	[det.n_moddefs] array with the number of cells in each module type
   det.cellcoords	[2 * max(n_cells) * det.n_moddefs] array with the cell (u,v) coordinates for each module type. The cell coordinates are defined relative to the module center, along the u and v coordinates, in mm.
   det.n_samples	[det.n_moddefs] vector with the number of samples per cell.
   det.samplecoords	[2 x max(det.n_samples) * det.n_moddefs] array of sample centers for each module type. All cells in a module type have the same sampling pattern. The sample coordinates are defined relative to the cell centers, along the u and v coordinates, and in mm. The samples may include x-ray crosstalk samples in neighboring cells.
   det.weights	[max(det.n_samples) * det.n_moddefs] array of relative sample weights for each module type. The weights have to sum up to one!
   det.activearea	[det.n_moddefs] array of the cell active area for a given module type.
   det.width	[det.n_moddefs] array of module widths. This serves as bounding box for all cell samples.
   det.height	[det.n_moddefs] array of module heights. This serves as bounding box for all cell samples.

   
   Detector definition

   Pack-Module structure and Multi-Z-Module
   

   Pack-Module structure and Multi-Z-Module are the new design of detector in Release 6.0. The detector is built as 2D arrays of modules, and a module is a 2D array of packs, a pack is a 2D array of cells, as shown in Figure 16, the detector is defined with 4 by 3 modules, the module is defined with 2 by 2 packs and the pack is defined with 3 by 3 cells.
   

   • Kerf size
     
   For each cell, the pitch size and kerf size are defined along both column and row direction. The kerf size is the width of the kerf between two cells, which is filled with reflector material, so that the active area (fill fraction) is defined.
   
   
   • Pack gap and module gap
     
   Pack gap/module gap are the gaps between two packs or two modules along u/v vector. This module allows defining different pack/module gap values for all the gaps.
   
   
   • Module rotation around u vector and position offset along w vector
     
     

   All the modules can have a rotation around u vector, and a position offset along w vector (w vector is the orthogonal vector of u and v vector). With that we can define the “chiclet” detector.
   

   
   Figure 16: Detector defined with Pack-Module structure

   
   

   The cfg variables are:
   

   % DEFINE PACK
   cfg.col_size = 1.0915;		                    % Detector column size or x-pitch (mm)
   cfg.row_size = 1.0915;			            % Detector row size or z-pitch (mm)
   cfg.col_cast = 0.120;			            % Reflector thickness at column direction (mm)
   cfg.row_cast = 0.120;			            % Reflector thickness at raw direction (mm)
   cfg.cols_per_pack = 8;			            % Number of detector columns per pack
   cfg.rows_per_pack = 16;			            % Number of detector rows per pack
   
   % DEFINE MODULE
   cfg.col_packs_per_mod = 2;		            % Number of column packs per module
   cfg.row_packs_per_mod = 2;		            % Number of row packs per module
   cfg.col_pack_gap = 0;			            % Inter-pack gap at column direction (mm)
   cfg.row_pack_gap = 0.12;		            % Inter-pack gap at row direction (mm) 
   
   % DEFINE DETECTOR
   cfg.mods_per_det_x = 52;                            % Number of modules per detector along X axis
   cfg.mods_per_det_z = 8;			            % Number of modules per detector along Z axis
   cfg.mod_gap_u = 0.;			            % Inter-module gap at column direction (mm)
   cfg.mod_gap_v = 0.1;			            % Inter-module gap at column direction (mm)
   cfg.col_offset = -4.25;			            % In-plane detector offset (in mm)
   cfg.mod_offset_w = zeros(1,cfg.mods_per_det_z);	    % module position offset along w vector
   cfg.mod_rotation_u = ones(1,cfg.mods_per_det_z)*0;  % module rotation around u vector
   cfg.mod_row_offset = zeros(1,cfg.mods_per_det_x);   % module position offset along row direction
   
   

   Sampling

   There are a number of variables that can be used to adjust sampling, in order to reduce simulation time by sacrificing precision, or increase precision at the expense of simulation time. These are:
   

   cfg.number_Ebins	number of energy bins used for the spectrum (1=monochromatic, empty=use the number of energy bins used in the source spectrum file)
   cfg.focal_oversample_x	number of samples in the x-direction for the focal spot.
   cfg.focal_oversample_y	number of samples in the y-direction for the focal spot.
   cgf.col_oversample	number of samples in the column-direction for the detector cells.
   cfg.row_oversample	number of samples in the row-direction for the detector cells.
   cfg.subviews_per_view	number of sub-samples for each view (at different rotation increments). cfg.subviews_per_view=1 represents a pulsed X-ray source.

   Note: In order to take full advantage of the focal spot and detector models, and correctly model spatial resolution, appropriate oversampling must be utilized. The required sampling depends on a number of factors, including focal spot specification (size and distribution), system geometry, phantom specification (feature size and location within the FOV), and detector specification (cell size and crosstalk). As a starting point, if simulating a 50 µm wire, the following oversampling is suggested.
   

   cfg.focal_oversample_x = 10;   % vary in the range of 7-15 to evaluate
                                  % dependency on this parameter.
   cfg.focal_oversample_y = 1;
   cgf.col_oversample = 30;       % vary in the range of 20-40 to evaluate
                                  % dependency on this parameter.
   cgf.col_oversample = 1;
   cfg.subviews_per_view =  1;    % vary in the range of
                                  %  1 for features at the center of the FOV,
                                  % 10 for features at the edge of the FOV.
   

   
   

   Detection Process

   The following variables specify the detector performance.
   

   Detection efficiency
   

   cfg.callback_detresponse	callback for detector behavior.
   cfg.detector_prefilter	series of materials and thicknesses to define pre-detector filtration.
   cfg.callback_dqe	callback for DQE.
   cfg.scintillator_material	name of the scintillator material. Default: cfg.scintillator_material = 'cesiumiodide';
   cfg.scintillator_depth	depth of the scintillator (mm). This impacts the detector efficiency. Use 1000 for 100% efficiency.

   Noise
   

   Energy weighting
   

   Lag
   

   cfg.callback_lag	computes lag. Default:cfg.callback_lag='Detection_Lag';
   cfg.lag_constant	fraction of each view that is added to the next view due to lag, afterglow, or other memory effects. Use 0.0 to eliminate lag.
   cfg.lag_compute	boolean that specifies whether or not to include lag in the simulations

   Optical crosstalk
   

   cfg.col_crosstalk_opt	fraction of the signal that is detected as optical cross-talk in neighboring detector columns. For example 0.05 means that each channel receives 5% of the x-rays from each of its neighboring channels. No longer used.
   cfg.row_crosstalk_opt	fraction of the signal that is detected as optical cross-talk in neighboring detector rows.

   DAS

   The following variables specify the DAS performance.
   

   cfg.callback_DAS	callback function for DAS.
   cfg.callback_DAS_air	callback function for DAS during air scans.
   cfg.callback_DAS_offset	callback function for DAS during offset scans.
   cfg.das_offset	offset for DAS.
   cfg.das_enoise	sigma of the electronic noise (in number of electrons) ; 0 means no electronic noise.
   cfg.das_gain	number of electrons collected for each 1 keV X-ray detected.
   cfg.das_lsb	least significant bit. The DAS will quantize the output based on this number. Use 0 to eliminate quantization.
   cfg.das_minvalue	Minimum value to truncate the DAS output. Use 0 to enforce positive numbers. Use for example -1e12 to avoid truncation.

   Reconstruction

   Third-generation clinical CT scanners employ a fan-beam or cone-beam geometry with a curved detector (figure 3, left-hand diagram), whereas projection x-ray systems (i.e. for radiography, diagnostic/interventional cardiology and neuroradiology, and mammography including 3D tomographic mammography), preclinical micro-CT and cone-beam CT (CBCT) scanners employ flat-panel detectors (figure 3, right-hand diagram). XCIST supports simulation of systems using either geometry, and reconstruction of projections from CT systems using either geometry. State-of-the-art CT scanners use multi-row detectors; the resulting cone-beam geometry requires 3D reconstruction.
   

   
   Figure 3: CatSim Detector-Based Coordinate Systems

   Tang’s 3D weighting algorithm [1] is an improved version of the classical Feldkamp-type cone-beam algorithm [2,3]. This approximate yet practical algorithm is well tested on 16- and 64-slice CT scanners, offering a balanced combination of low cone-beam artifacts with good noise and spatial resolution performance. While its general formulation applies to helical cone-beam CT, it extrapolates well to the case of axial cone-beam scanning. Axial FDK-based algorithms are very popular and widely applied in clinical and preclinical applications. We include a few sample reconstruction kernels including ‘soft’, ‘standard’, and ‘bone’, which represent typical clinical reconstruction kernels. The modular design enables users to use custom kernels. Parker weighting [4] is used for short-scan simulations

   UML Reconstruction

   Prerequisites: Install XCIST using pip install gecatsim.

   How the recon code works: The XCIST reconstruction code utilizes the FDK-type algorithm and supports both axial and helical scan modes. Implemented in C/C++ for enhanced performance, it can be accessed via a Python interface. In axial reconstruction, a ramp filter is applied first, followed by back projection. For helical reconstruction, the process involves fan-to-parallel rebinning of the sinogram data, followed by ramp filtering and then back projection. For further details, please refer to references shown below [5,6].

   Tunable Parameters: The following parameters can be specified in XCIST to control the reconstruction process:
   

   recon.fov                               # diameter of the reconstruction field-of-view (in mm)
   recon.imageSize                         # number of column/row pixels to be reconstructed (currently only square is supported)
   recon.sliceCount                        # number of slices to reconstruct
   recon.sliceThickness                    # reconstruction slice thickness AND inter-slice interval (in mm)
   recon.centerOffset = [0.0, 0.0, 0.0]    # reconstruction offset relative to center of rotation (in mm)
   recon.reconType = 'fdk_equiAngle'       # Name of the recon function to call, can be fdk_equiAngle for axial scan or helical_equiAngle for helical scan
   recon.kernelType = 'standard'           # reconstruction kernels. 'R-L' for the Ramachandran-Lakshminarayanan (R-L) filter, rectangular window function,'S-L' for the Shepp-Logan (S-L) filter, sinc window function, 'soft', 'standard', 'bone' for kernels similar to those on clinical scanners
   recon.startAngle = 0                    # in degrees; 0 is with the X-ray source at the top
   recon.unit = 'HU'                       # unit of recon images, can be '/mm', '/cm', or 'HU'
   recon.mu                                # mu in /mm to convert recon to HU; typically around 0.02/mm
   recon.huOffset                          # unit is HU, -1000 HU by definition but sometimes something else is preferable
   recon.printReconParameters = False      # Flag to print the recon parameters
   recon.saveImageVolume = True            # Flag to save all recon slices as one big file
   recon.saveSingleImages = False          # Flag to save recon results as individual images
   recon.printReconParameters = False      # Flag to print the recon parameters
   recon.displayImagePictures = False      # Flag to display the recon results as .png images
   recon.saveImagePictureFiles = False     # Flag to save the recon results as .png images
   recon.displayImagePictureAxes = False   # Flag to display the axes on the .png images
   recon.displayImagePictureTitles = False # Flag to display the titles on the .png images
   

   How to run recon in python: Here is an example Python script to run the reconstruction, assuming the desired sinogram file “{ct.resultxName}.prep” is already in the same folder. The “{ct.resultxName}.prep” file should contain binary float32 data in little endian format with dimensions [number of views, number of detector rows, number of detector columns]. Example configuration files can be downloaded from: https://github.com/xcist/example/tree/main/AAPM_datachallenge/simulation_scripts.
   

   import gecatsim as xc
   from gecatsim.reconstruction.pyfiles import recon
   # import cfg files
   ct = xc.CatSim("AAPM_MAR_phantom",
                            "AAPM_MAR_protocol",
                            "AAPM_MAR_physics",
                            "AAPM_MAR_scanner",
                            "AAPM_MAR_recon")
   ct.resultsName = "out"
   
   # reconstruction
   cfg = ct.get_current_cfg();
   # so xcist knows we are doing recon
   cfg.do_Recon = 1
   
   # run recon
   recon.recon(cfg)
   
   

   The following variables specify the reconstruction process.
   

   cfg.recon_fov	reconstruction field of view.
   cfg.recon_size	size of reconstructed image (default = 512x512).
   cfg.recon_planes	number of “slices” (XY planes) in the z dimension.
   cfg.recon_slice_thickness	thickness of a slice during reconstruction.
   cfg.recon_xcenter	center of reconstruction in x dimension.
   cfg.recon_ycenter	center of reconstruction in y dimension.
   cfg.recon_zcenter	center of reconstruction in z dimension.
   cfg.recon_filter	filter used during reconstruction.
   cfg.recon_type	type of reconstruction algorithm.
   cfg.recon_mu	linear attenuation coefficient used in the conversion to Hounsfield units.
   cfg.recon_unit	specifies the unit in which reconstructed images are output (HU, 1/cm).
   cfg.recon_hu_offset	specifies the offset applied after normalization (typically 0 or -1000).
   cfg.recon_fbp3d_write_intermediate_results	1 = save a slice after backprojecting each view, & write to file.
   cfg.recon_fbp3d_intermediate_result_slice_index	slice for which to save intermediate results (0 = middle slice).

   Gantry geometry

   The overall gantry geometry is by
   

   cfg.sid = Source-To-IsocenterDistance_Millimeters;
   cfg.sdd = Source-To-DetectorDistance_Millimeters;
   

   The detector can be positioned slightly off-center relative to isocenter. This is specified by
   

   cfg.col_offset = ColumnOffset_Columns;
   cfg.row_offset = RowOffset_Rows;
   

   which give the relative offset in the u and v direction respectively, in units of columns and rows.
   

   This can be used to achieve quarter detector offset.
   

   The lateral center of the detector, in units of columns, is

   cfg.col_center = (1 + cfg.col_count)/2.0 - cfg.col_offset;
   

   The longitudinal center of the detector, in units of rows, is

   cfg.row_center = (1 + cfg.row_count)/2.0 - cfg.row_offset;
   

   Scan protocol

   The gantry rotation and table translation are defined by the gantry callback function:
   cfg.callback_gantry = 'Gantry_Helical';
   The same function can be used to compute vibration, multi-source acquisitions, focal spot wobble, etc. The number of views per 360 degrees of rotation is defined by
   cfg.views_per_rotation = NumberOfViews;
   

   The tube current (without mA modulation) is given by
   cfg.mA = XrayTubeCUrrent_mA;

   The X-ray quantum noise can be disabled by setting
   cfg.disable_xray_noise = XrayNoiseDisableFlag;

   The seed of the random generators can be set by
   cfg.seed = []; or
   cfg.seed = SeedInteger;
   If the seed is empty, every run has a different noise realization. If it is set to a particular (integer) value, the same random realization will be used every time.
   

   The mA modulation callback routine and level are set using
   cfg.callback_mA_modulation = []; or
   cfg.callback_mA_modulation = 'mA_modulation_sin';
cfg.mA_modulation = 0;
   If the callback function field is empty, no mA modulation is applied. Otherwise, a level of 0.0 means NO modulation and 1.0 means modulation between 0% and 200%.
   

   The gantry rotation speed is defined by
   cfg.rotation_period = GantryRotationPeriod_Seconds
   

   The initial time (in s, usually 0s), the initial table position (in mm) and the initial rotation angle (in degrees) of the simulation are defined by
   cfg.start_time = ScanStartTime;
cfg.start_z = ScanStartZ_Millimeters;
cfg.start_angle = ScanStartAngle_Degrees;
   

   For example, if an object is positioned at -100mm, it will be at the center of the gantry for a table z-position of +100mm. The table speed (in mm/s) is defined along the positive z-axis by
   cfg.table_speed = TableSpeed_MillimetersPerSecond;
   which is zero for axial acquisitions.
   

   The rotation direction is defined by
   cfg.rotation_direction = GantryRotationDirecton;
   where 1 means counterclockwise and –1 means clockwise (as viewed from the back) of the gantry.
   

   The total scan length can be defined by the total number of views (always integer) or by the total number of rotations (can be non-integer):
   cfg.total_n_views = 5*cfg.views_per_rotation;
   or
   cfg.n_rotations = 5;
   

   The gantry tilt angle (in degrees) is defined by
   cfg.tilt_angle = GantryTiltAngle_Degrees;
   The gantry is tilted around the x-axis: positive tilt angle means that for positive y, the mid-plane is tilted away in the negative z-direction (left-handed convention).
   

   It is possible to break up simulations into a number of sub-simulations by defining particular start and stop views:
   cfg.start_view = StartViewNumber;
cfg.stop_view = StopViewNumber;
   Normally a simulation would start at view 1 and stop at view cfg.total_n_views.
   

   Air and offset scans typically consist of a single noise-free view (in simulations):
   cfg.airscan_total_n_views = AirScan_NumberOfViews;
cfg.offsetscan_total_n_views = OffsetScan_NumberOfViews;
   except for multi-source simulations in which each source produces a different air flux at the detector.
   

   Afterglow

   Prep

   The following variables specify the prep process.
   

   cfg.callback_pre_log	specifies the names of the post-processing routines to be applied before the log-conversion step
   cfg.callback_log_conversion	specifies the routine to do the log-conversion
   cfg.callback_post_log	specifies the names of the post-processing routines to be applied after the log-conversion step

   Scatter (convolution model)

   The routine to add scatter to the simulation is defined by
   

   cfg.callback_scatter = 'Scatter_ConvolutionModel';
   

   If cfg.callback_scatter is empty ([]), no scatter will be simulated. ScatterConvolutionModel uses a simple sinogram-based scatter estimation and requires the following parameters :

   cfg.scatter_level = 0.001;
   cfg.scatter_smoothing = 10;
   

   cfg.scatter_level is the scatter fraction relative to the airscan, i.e., scatter-to-air ratio (SAR), NOT scatter-to-primary ratio (SPR). Using 0.0 will eliminate scatter.
   

   SAR should be on the order of 0.001 or so. The reason it was not defined by SPR is that it was not completely obvious how to define the primary level.
   

   cfg.grid_height specifies the height of anti-scatter grid. This is only used for Monte Carlo reconstruction.

   Dose Reconstruction

   Dose Reconstruction is a new method for volumetrically reconstructing absorbed dose on a per-voxel basis, directly from the actual CT images.
   

   The interface of Dose Recon function is very similar as CatDose.
   

   The main function is catdoserecon.m:
   

   function [dosevol] = catdoserecon(configfilename, cfg, adjust, preadjust, silent)
   

   The cfg variables are:

   cfg.doserecon_n_voxel = cfg.recon_size;		  % size of the volume for dose recon
   cfg.doserecon_xoffset = 0;			  % patient x offset from iso-center
   cfg.doserecon_yoffset = 0; 			  % patient x offset from iso-center
   cfg.doserecon_zoffset = 0; 			  % patient x offset from iso-center
   cfg.doserecon_Mu_water = 0.2; 			  % effective Mu according to spectral cal
   cfg.doserecon_Mu_bone = 0.75; 			  % effective Mu according to spectral cal
   cfg.doserecon_do_convol = 0; 			  % do convolution or not (1/0)
   cfg.doserecon_phantom_material = 'polyethylene';  % if phantom material is known ...
   

   • How to use:
     

        I. Generate the CT image.
        II. Call catdoserecon function.
   

   For example, with catsim and catsimrecon:

   catsim xxx.cfg;
   catsimrecon xxx.cfg;
   catdoserecon xxx.cfg;
   

   References:
   

   [1] Tang X et al 2006 "A three-dimensional-weighted cone beam filtered backprojection (CB-FBP) algorithm for image reconstruction in volumetric CT-helical scanning." Phys. Med. Biol. 51 855–74
   [2] Feldkamp L A, Davis L C and Kress J W 1984 "Practical cone-beam algorithm J. Opt. Soc. Am." A 1 612
   [3] Wang G, Ye Y and Yu H 2007 "Approximate and exact cone-beam reconstruction with standard and non-standard spiral scanning." Phys. Med. Biol. 52 R1–13
   [4] Parker D 1982 "Optimal short scan convolution reconstruction for fanbeam CT" Med. Phys. 9 254–7
   [5] Wu, Mingye, et al. "Xcist—an open access x-ray/ct simulation toolkit." Physics in Medicine & Biology 67.19 (2022): 194002.
   [6] Tang, Xiangyang, et al. "A three-dimensional-weighted cone beam filtered backprojection (CB-FBP) algorithm for image reconstruction in volumetric CT—helical scanning." Physics in Medicine & Biology 51.4 (2006): 855.