31. X-ray Imaging Unit

The X-ray Imaging unit obtains an X-ray image of the domain at specific simulation time(s) using the main driver function XrayImaging. The implementation is based on sending one X-ray per detector screen resolution pixel through the domain and record it's intensity reduction due to attenuation in the domain media. Currently only cold atomic opacities are being used to determine the combined opacity of each cell's materials. Each detector can record X-ray intensities for a mix of X-ray energies, thus each X-ray for each detector in the code represents a bundle of X-rays with predefined distributions of X-ray energies. Each detector is also given an X-ray point origin from which the X-rays originate.

31.0.1 X-ray Paths through Domain

All X-rays travel in straight lines through the domain. Each X-ray is given a dimensionless 3D unit direction vector ${\bf u},\vert{\bf u}\vert=1$ with components $(u_x,u_y,u_z)$ pointing from the detector screen to its X-ray origin. Each direction vector component is a direct measure of how much 3D distance the X-ray travels between two such component points. For example, the direction vector component $u_x$ can be used to evaluate the 3D distance $d$ covered by the X-ray between two points ${\bf p}_1$ and ${\bf p}_2$ by using their components $p_{1x}$ and $p_{2x}$:

$$d = \vert(p_{1x}-p_{2x})/u_x\vert,$$ (31.1)

where, when removing the absolute value bars, we can judge, if the X-ray actually travelled that distance (positive $d$) or not (negative $d$), depending on the location of $p_{1x}$ and $p_{2x}$ and the sign of $u_x$. For 3D cartesian domains it is therefore straighforward to find the minimum distance travelled through a cell by chosing the minimum of all positive distances that are obtained from the above Eq.(31.1) when replacing $p_{1x}$ and $p_{2x}$ by the 3D cartesian cell bounding box component pairs. Once $d$ is found, we can use it to find the new X-ray's location ${\bf p}_{new}$ on the cell's bounding box from its old location ${\bf p}_{old}$:

$${\bf p}_{new} = {\bf p}_{old} + d * {\bf u}.$$ (31.2)

This simple procedure for 3D cartesian domian geometries changes for 3D X-ray tracing in 2D cylindrical domain geometries. While the axial z-component is identical to the 3D cartesian case, the curved mantle in the radial space makes it necessary to solve quadratic equations. Eq.(31.2) can be seen as the parametric line equation of a 3D X-ray. We only look at the x- and y-components of the line equation and add their squares together:

$$r_{new}^2 = r_{old}^2 + 2(p_{old,x}u_x+p_{old,y}u_y)d + (u_x^2+u_y^2)d^2,$$ (31.3)

where

$$r_{old}^2 = p_{old,x}^2+p_{old,y}^2$$ (31.4)

$$r_{new}^2 = p_{new,x}^2+p_{new,y}^2$$ (31.5)

are the old and new radial 2D cylindrical positions of the 3D X-ray. The quadratic equation in $d$ is solved for positive $d$ by inserting the 2D cylindrical cell radial bounding box pair $(r_{min},r_{max})$. Due to the necessity of solving quadratics, the computational accuracy attainable for the radial $d$ can be lower than for the axial $d$, which is solved by using the linear Eq.(31.1). Instead of forcing the quadratic $d$ solution to be as accurate as the linear (axial) $d$ solution (for example by solving the quadratic in higher precision), the linear X-ray propagation through the cell is allowed to take place using more than one step (at most two steps are ever required). The 3D X-ray tracing in a 2D cylindrical domain also requires constant adjustments between its 3D x- and y-position components and its 2D radial component, as for example cell nudging of the X-ray must occur in the 2D cylindrical radial space. One cannot allow the X-ray's x- and y-position components to get out of sync with the radial component, as all of them are needed to solve Eq.(31.3).

31.0.2 Creating and Tracing the X-rays

Each pixel center coordinate on the detector screen is the starting point of an X-ray. It's trajectory is a straight line from the pixel center to the X-ray point origin of the detector screen. The X-ray's direction vector ${\bf u}$ is pointing from the pixel center to the origin. Initial domain entry points for each X-ray are determined using their ${\bf u}$ and Eq.(31.2), where ${\bf p}_{old}$ is the X-ray's global pixel center coordinates and ${\bf p}_{new}$ the domain boundaries. If an X-ray misses the domain, it is recorded directly as a screen X-ray with zero attenuation. There is no screen miss possibility for the X-ray, because per construction each X-ray belongs to a pixel on the detector screen. If the X-ray hits the domain, its entry position and direction vector are stored to prepare for X-ray tracing through the domain, as well as it's initital attenuation value is set to zero. Besides position, direction and attenuation information, each domain X-ray needs additional data to help with its movement through the domain and its transformation to a screen X-ray, such as current block and processor numbers, a global tag number and the integer pixel pair label from where it originated on the detector screen. The screen X-rays will only possess the local screen coordinate pair, it's calculated relative intensity in the range $[0,1]$ and the detector screen number.

The mechanism for tracing the X-rays through the domain resembles the one used for tracing protons through the domain. X-rays are created and traced through a frozen domain at a particular time step, hence no time resolved X-ray tracing is needed. Batches of domain X-rays are generated and emptied into screen X-ray buckets once the X-rays leave the domain. An important difference when compared to proton imaging is that each batch of X-rays must correspond to one specific detector screen. The reason behind this is that calculation of the cell opacities depends on the X-ray energies used and have to be calculated for each detector separately. The X-ray energies hence influence the domain properties through material opacities, whereas the proton energies only have an influence on their initial velocities and this property will not change any domain values.

31.0.3 Setting up and Recording on the X-ray Detector Screens

The X-ray detector screens are identical in shape and specifications to the ones used for proton imaging (see section 28.0.4). The same parameters are used to specify size (sidelength), location (global position ${\bf C}$ of the detector center) and orientation (screen normal vector ${\bf n}$ and tilting axis and angle) in space. Additional parameters used specifically for X-ray detectors are: 1) global origin ${\bf O}$ position as the source point for the X-rays, 2) activation time, indicating the simulation time when the X-rays should be recorded, 3) screen resolution $R$, which indicates how many pixels will be along each detector side (leading to a total of $R^2$ pixels), 4) the energy distribution for each X-ray (finite set of X-ray energies and fractions of these X-ray energies adding up to 1) and 5) data related to each pixel size to facilitate recording of the X-rays on the screen. Recording the X-rays on the detector screen is much simpler than for protons, since the target pixel is already known. The X-rays integer pixel coordinate pair is converted to the local screen coordinate pair and the accumulated attenuation $a$ is translated to it's relative screen intensity $I_{rel}$:

$$I_{rel} = I_{exit}/I_0 = I_0e^{-a}/I_0 = e^{-a},$$ (31.6)

where $I_{exit}$ is the X-ray intensity when exiting the domain and $I_0$ the initial X-ray intensity. Note, that currently the initial intensity needs not be supplied as part of the detector info, since only relative intensities are recorded. As for proton imaging, each X-ray detector file is set up as a formatted ascii file with name:

$$\mbox{$<$basename$>$XrayDetectorFile$<$detectorID$>$\_$<$timestamp$>$}$$ (31.7)

where $<$basename$>$ is the simulation name, $<$detectorID$>$ is the detector number (currently limited to the range [01-99]) and (optionally) $<$timestamp$>$ is the simulation time when the file was created. Three quantities are written to the detector file for each X-ray: 1) the local screen coordinate pair $(x,y)$ in the range $[0,1]$ and 2) it's relative intensity $I_{rel}$. Translation to a greyscale PGM picture can be done by using the supplied ASCII $->$ PGM greyscale converter software.

31.0.4 Usage

The following additions to the setup command activate the X-ray Imaging unit:  

+XrayImaging [xray_maxDetectors=<number> threadXrayTrace=True]

The +XrayImaging shortcut handles all the logistics for properly including the XrayImaging unit. The setting options are:

• xray_maxDetectors: The maximum number of X-ray imaging detectors for the simulation (default = 1).
• xray_maxEnergyLevels: The maximum number of X-ray energy levels for a detector (default = 3).
• threadXrayTrace: Enables threading during X-ray domain tracing (default = false).

The following is a list of flash.par runtime parameters for the X-ray imaging unit.

31.0.4.1 X-ray Imaging Detectors Runtime Parameters

The following are the runtime parameters for the X-ray imaging detectors. The _n at the end of each runtime parameter characterizes the detector number and needs to be replaced by _1 for the first detector, _n by _2 for the second detector, etc.

• xray_numberOfDetectors: The number of X-ray imaging detectors that are going to be used.
• xray_detectorCenter[X,Y,Z]_n: The global 3D components of the detector center position vector ${\bf C}$.
• xray_detectorOrigin[X,Y,Z]_n: The global 3D components of the X-rays origin position vector ${\bf O}$.
• xray_detectorNormal[X,Y,Z]_n: The local 3D components of the detector normal vector ${\bf n}$.
• xray_detectorEnergyLevelCount_n: The number of X-ray energy levels.
• xray_detectorEnergy_n_Level_i: The i-th X-ray energy level in eV. There must be as many X-ray energy levels as the number of X-ray energy levels specified.
• xray_detectorEnergy_n_Fraction_i: The i-th X-ray energy level fraction contribution. There must be as many X-ray energy level fractions as the number of X-ray energy levels specified. All fractions must be in the range [0-1], and they must add up to 1.
• xray_detectorPerpXrays_n: If this parameter is set true, the origin X-ray position will be ignored and all the X-rays will be hitting the screen in a perpendicular fashion, i.e. parallel to the screen's normal direction vector. This is to avoid having to specify origin positions mimicking an infinite distance away from the screen.
• xray_detectorResolution_n: The number of pixels along each side of the detector.
• xray_detectorSideLength_n: The side length (cm) of the detector square screen.
• xray_detectorSideTiltingAngle_n: The tilting angle (degrees) of two sides parallel to the detector screen unit axis ${\bf u}_y]$ with respect to the tilting axis.
• xray_detectorSideTiltingAxis_n: The global tilting axis ('x','y' or 'z') for the detector screen.
• xray_detectorAlignWRTorigin_n: A useful shortcut to orient the detector screen in such a way that the normal vector coincides with the origin / detector center line. Given the origin ${\bf O}$ and the detector center ${\bf C}$, the normal vector ${\bf n}$ is recalculated such that $({\bf O}-{\bf C})\cdot {\bf n} = 0$.

31.0.4.2 X-ray Imaging General Runtime Parameters

• xray_cellWallThicknessFactor: Controls the (imaginary) thickness of the cell walls to ensure computational stability of the X-ray imaging code. The cell thickness is defined as this factor times the smallest cell dimension along all geometrical axes. The factor is currently set to $10^{-6}$ and should only very rarely be changed.
• xray_3Din2D: Indicates, whether this is a 3D X-ray in 2D cylindrical domain simulation.
• xray_detectorFileNameTimeStamp: If set to true, a time stamp will be added to each detector file name. This allows for time splitting of detectors.
• xray_detectorXYwriteFormat: Controls formatted ascii output of the screen X-ray local coordinate $(x,y)$ pairs on the detector screens (default 'es20.10').
• xray_detectorDGwriteFormat: Controls formatted ascii output of the screen X-ray relative intensities $I_{rel}$ on the detector screens (default 'es15.5').
• xray_maxXrayCount: The maximum number of X-rays that can be created on one processor per batch in the domain.
• xray_printDetectors: If true, it prints detailed information about the X-ray detectors to a file with name $<$basename$>$XrayImagingDetectors.txt, where $<$basename$>$ is the base name of the simulation.
• xray_printMain: If true, it prints general information regarding the X-ray imaging setup to a file with name $<$basename$>$XrayImagingMainPrint.txt, where $<$basename$>$ is the base name of the simulation.
• xray_printSpecies: If true, it prints general information regarding the all the species present in the current simulation to a file named $<$basename$>$XrayImagingMainSpecies.txt, where $<$basename$>$ is the base name of the simulation. It records the atomic element decomposition of each species with number and mass fractions and also records the cold opacities of each species for each detector (depending on the X-ray energy specified).
• xray_printXrays: If true, it activates printout of domain X-ray data at specific stages during the domain X-ray tracing. This runtime parameter is meant solely for debugging purposes and calls to the corresponding X-ray printing routine can be placed anywhere in the X-ray imaging code. Each processor writes its own file(s) with name(s):
  $<$basename$>$printXraysDetBatch$<$DetectorBatchNr$>$Proc$<$ProcessorID$>$_$<$simulationTime$>$.txt,
  where $<$basename$>$ is the base name of the simulation. An optional _Tag$<$Tagnumber$>$ can be added to these files to distinguish between printouts at different strategic places. Use of this feature is reserved ONLY for debugging purposes and is currently limited to 1000 detector batches and 1000 processors per time stamp. Users other than code developers should not activate this feature.
• xray_XrayDeterminism: If true, the Grid Unit will be forced to use the sieve algorithm to move the X-ray particle data. Forcing this algorithm will result in a slower movement of data, but will fix the order the processors pass data and eliminate round off differences in consecutive runs.
• xray_screenXrayBucketSize: Sets the bucket size for flushing screen X-rays out to disk.
• useXrayImaging: If false, no X-ray imaging will be performed, even if the code was compiled to do so. Bypasses the need to rebuild the code.
• threadXrayTrace: If true, X-ray tracing through a block is threaded. This runtime parameter can only be set during setup of the code.

31.0.5 Unit Test

The unit test supplied for the X-ray imaging unit is based on putting a sphere of pure iron inside a domain filled with pure hydrogen. The hydrogen opacity is enforced to be equal to zero, hence attenuation of the X-rays is solely due to the atomic opacity of iron in the sphere. For each type of domain (3D cartesian and 2D cylindrical), three detector screens, each with resolution R = 1000, are placed at three different positions in space. For the cubic 3D cartesian domain, the three screen positions are along a face center, edge center and vertex direction. For the 2D cylindrical domain, we also place the three screens along the radial, axial and cylinder edge of the implied 3D cylinder. Attenuation of the X-rays are recorded and compared to analytically obtained bounds.

31.0.5.1 X-ray Attenuation inside a Sphere

Given a sphere with radius $R$ and located with its center at ${\bf C}$, any X-ray with screen location ${\bf S}$ and origin location ${\bf O}$ will travel a distance $d$ inside the sphere given by:

$$d = 2\sqrt{\max(0,R^2-\sin^2\phi * \vert CS\vert^2)},$$ (31.8)

where $\vert CS\vert$ is the distance between ${\bf C}$ and ${\bf S}$ and $\phi$ is the angle between the vectors $\overrightarrow{\bf SC}$ and $\overrightarrow{\bf SO}$. The sphere representation inside the domain is blocky in all three directions for 3D cartesian domains and blocky in the axial direction (stacked cylindrical disks) for 2D cylindrical domains. Hence Eq.(31.8) cannot be directly used to determine the exact distance travelled through the blocky spheres. However, one can define two extra spheres: 1) Maximum inner sphere, with radius $R_{min}$, which is the largest sphere that can be placed completely inside the blocky spheres and 2) Minimum outer sphere, with radius $R_{max}$, which is the smallest sphere that can contains the blocky spheres completely. Since the positions of the X-ray don't change, $\phi$ and $\vert CS\vert$ does not change and each $R_{min}$ and $R_{max}$ also define through Eq.(31.8) corresponding X-ray distances $d_{min}$ and $d_{max}$ travelled through the maximum inner and minimum outer spheres. We then have the following three inequalities

$$R_{min}\leq R \leq R_{max}$$

$$d_{min}\leq d \leq d_{max}$$ (31.9)

$$I_{min}\geq I \geq I_{max},$$

where $R$,$d$ and X-ray intensity $I$ correspond to the blocky sphere quantities. Using $d_{min}$ and $d_{max}$, upper and lower bounds for the X-ray attenuation can be given when crossing the domain, leading to X-ray intensities $I_{min}$ and $I_{max }$ due to the maximum inner and minimum outer spheres.

31.0.5.2 FLASH Setup of the Unit Test

For the 3D cartesian unit test, we set up a cubic box with side length 10 cm. A pure iron sphere with density 0.3 g/cm$^3$ and radius 5 cm is created at the center of the domain, using uniform refinement level 6 for the AMR grid. Detector screens with resolution R = 1000 and of side length 25 cm are placed at the following 3 different places, such that the sphere center and the detector screen centers are 10 cm apart: 1) along the y face center direction (0,1,0), along the edge center direction (1,1,0) and along the vertex direction (1,1,1). The corresponding X-ray origins are all placed in the opposite directions with origin / screen center distance of 20 cm. All X-ray energies are set to 60 keV. The three intensities $I$, $I_{min}$, $I_{max }$ of each X-ray from all three detectors are analyzed and the unit test is marked successful, if for all X-rays the intensity inequality in Eq.(31.9) is fulfilled.

The 2D cylindrical unit test is set up to match the 3D cartesian one. The 2D cylindrical domain has an axial dimension of 10 cm and a radius dimension of 5 cm. The implicit 3D sphere is placed in form of a semicircle with center at the domain location (radial,axial) = (0,5). All other sphere specifications (material and density) and detector specifications (side length 25 cm, R = 1000, detector center to detector origin distance = 20 cm, detector center to sphere center distance = 10cm) match the 3D cartesian unit test. The three detectors are placed as follows: 1) along the axial direction, 2) along the radial direction and 3) along the cylindrical edge direction. The 2D cylindrical unit test is successful, if all X-ray intensities are within the bounds mentioned previously.