This post will review the concepts behind medical imaging and the forward model used in X-ray imaging. Additionally, this tutorial will explain the development of the forward model in Python.
X-rays are one type of electromagnetic radiation used in medical imaging. The setup of an X-ray imaging systems consists of an X-ray transmitter (X-ray Tube) and a detector (Flat Panel Detector). X-rays interact with any matter between the transmitter and detector, for medical imaging this could be a femur, stomach or another area on the patient. As a result of an object in the imaging region, the transmitted X-rays lose energy due to photon interactions with matter and produce a projected pattern of this matter on the detector. One way to describe the amount of energy that is lost per unit length by the transmitted X-ray is which is the linear attenuation coefficient. The linear attenuation coefficient is different in every part of the body due to the different matter in each of these regions.
X-ray Imaging Model
One way to describe X-ray imaging is it produces a projection image of the combined linear attenuation coefficient onto the detector. This can be seen by the X-ray modeling equation of
This first term in the equation is which is the energy emitted from the X-ray source. The next term is the previously stated linear attenuation coefficient in a certain X, Y, Z coordinate and L is the path length. Eq. 1 is the forward model which describes the energy received at a detector pixel given an initial X-ray energy and discretized model for the linear attenuation coefficients. This forward model is important for testing out different X-ray system parameters such as beam energy, beam filtering, and many others.
Discretization is a process used to take a continuous space and break it down into a manageable subspace that can fit into the memory on a computer. The reason why this is necessary for our forward model is to make the problem computationally tractable and sampling for linear attenuation coefficient matrices are discrete.
One example of discretization can be seen in the Fig. 3, where we break up the imaging space into discrete cubes (seen in purple). Each of these chunks is an area of the body in the 2D case that represents the linear attenuation coefficient in that area. One factor to consider when performing discretization is the resolution of each of these chunks. In the below figure, some chunks have both bone and air present which will result in a fuzzy representation of this area. When performing imaging it is important to have the right discretization to best resolve the various areas of the patient.
The model that relates the energy between an X-ray source and a detector pixel is a recursive implementation of Eq. 1. To apply this model, ray tracing is performed where a ray is stepped from source to detector pixel(Fig. 3). This ray will move along its direction vector and when it hits a new discretization in the imaging region will lose energy according to Eq. 1. A factor into the energy lost is the path length which is the distance the X-ray traveled within the discretization. This factor comes from the linear attenuation coefficient which is energy loss per unit length, so this path length is used to describe the amount of energy lost in that chunk.
The below shows an example of the 3D model that could be used for the above figure. Initially, we have the main loop function which will calculate the energy at each pixel location. To perform this calculation first, the direction vector between a pixel and the source is calculated. A ray is then initialized at the source and moves towards the next discretization cube on the direction vector. The move cube function is what performs the movement from one cube to another and returns the position and distance to the next cube. For each of these moves, the position is verified to be in the imaging region if, so it will perform the energy loss calculation, otherwise, it will continue the chunk movement (since outside of the imaging region is air where no energy is lost). To perform this energy calculation, Eq. 1 is used.
The cube marching function will take the current position and find the next cube for which the ray will land. The first line in this function will calculate the distance to hit the different dimensions of the cube (X, Y, Z). Then, we want to use the minimum times as this is where the ray intersects first. After that, we want to record the distance to get to that new location and update the new position. Note this new position has an added epsilon for the next update because some next iterations result in endless looping.
After performing all the above steps, one example of using these functions on a skull’s linear attenuation coefficients can be seen below:
For full code/constants definitions/imaging model diagram see Github, this code also has optimizations present which you can read about here.
That wraps up the basic model of X-ray imaging. One problem with the above Python implementation is it takes a long time to perform the forward model. In my next blog post, I will review Numba a library to accelerate your Python code and apply it to the X-ray imaging problem.