Iterative Image Reconstruction Methodology in Optical CT Radiochromic Gel Dosimetry

Modern advancements in radiation therapy require paralleled advancements in the dosimetric tools used to verify dose distributions. Optical computed tomography (CT) imaged radiochromic gel dosimeters provide comprehensive, tissue equivalent, 3D dosimetric information with high spatial resolution and low imaging times. Traditional CT image reconstruction methods (filtered backprojection) do not account for light refraction within the optical CT system reducing the image quality. Iterative reconstruction methods make use of a system matrix that describes this light refraction thus, improving the reconstructed image quality. However, use of iterative reconstruction methods is not widespread, largely due to the impractical storage size of the required system matrix. Furthermore, current iterative reconstruction methods do not address the issue of image degradation due to a single detector element collecting light from multiple raypaths. For optical CT radiochromic gel dosimetry to be used effectively as a radiation therapy treatment plan verification tool, the system must be both practical and accurate. Thus, this work has two main objectives: (i) reduce the size of system matrices by means of polar coordinate discretization in lieu of the traditional Cartesian coordinate discretization, and (ii) reduce image degradation due to multiple raypaths by a novel approach to populating the system matrix that accounts for multiple raypaths.


Introduction
Radiation therapy treatments plans with complex 3D dose distributions and steep dose gradients need dosimetric tools that can verify these complex treatment plans [1,2].Traditional dosimeters are restricted to point or planar measurements, whereas, gel dosimeters can measure dose in full 3D geometry while being nearly tissue equivalent [3,4].Furthermore, their radiation induced response can be measured with sub-millimetre spatial resolution [1,5].Specifically, radiochromic gel dosimeters elicit an optical density change proportional to the delivered radiation dose, allowing for readout of the dose distribution via optical computed tomography (CT).For optical CT radiochromic gel dosimetry to be used in radiation therapy treatment plan verification the system needs to be both accurate and practical.That is, the optical CT images must give an accurate representation of the dose distribution as well as being practically viable for a clinical setting.

System Matrix
Light is refracted as it traverses through different mediums in the optical CT system causing image degradation and artifacts when reconstructed using traditional filtered backprojection (FBP) methods.FBP assumes a straight raypath from source to detector, which is not the case when light undergoes refraction through the optical CT system.This loss of image quality can be remedied by utilizing iterative image reconstruction methods that account for this refraction.The path of each light ray as it travels through the optical system is compiled into what is known as a system matrix and is applied during iterative image reconstruction.These system matrices can become immensely large (>100GB) making them impractical for storage on a typical clinical workstation computer [6,7].Here we develop a method to reduce the storage size of optical CT system matrices.

Single Detector Collection of Multiple Raypaths
As light is refracted through an optical CT system there arises the potential for a single detector to collect light from multiple (often very different) raypaths.This violates the single raypath assumption of the linear attenuation equation, which is utilized in the iterative reconstruction, causing artifacts in the resulting image.Engineering efforts can be made during the design of an optical CT scanner to minimize this occurrence by means of lenses and refractive index match materials.Nonetheless, the problem of multiple ray paths persists as all interfaces, including the dosimeter material, would be required to be perfectly index matched which is difficult to achieve in practice.Some rays may undergo severe refraction such that they are collected by detectors that also received a ray from a more direct raypath.These severely refracted rays are termed crossover rays and the rays that follow a more direct raypath are termed primary rays.Furthermore, light is partially reflected at each medium interface which allows for the potential that a single ray can undergo multiple reflections before being collected by the detectors and are termed reflected rays.Regardless of whether a detector element receives multiple raypaths due to crossover or reflected rays the problem remains the same; single detector collection of multiple raypaths.Here we aim to address the issue of multiple raypaths entering a single detector element and thereby reduce the induced image artifacts.

Materials and Methods
An in-house ray tracing simulator has been developed to track the path of each ray as it traverses the different mediums of the prototype tabletop solid-tank fan-beam optical CT scanner developed at the University of British Columbia -Okanagan (Figure1) [8].The scanner consists of an acrylic (PMMA) block with a curved face and a hole bored through it in which a cylindrical radiochromic gel dosimeter is inserted.The dosimeter is submerged in a refractive index matched fluid to allow for the dosimeter to be rotated to acquire projections while the laser source and photodiode detector bay remain fixed.By discretizing the dosimeter into cartesian pixels, the simulation allows for the calculation of the length of every ray through every pixel and these calculated lengths were used to populate a system matrix.Furthermore, the ray tracing simulation was used to generate synthetic detector data (sinograms) mimicking the data collection of the actual optical CT system.These synthetic sinograms were used to evaluate the novel iterative reconstruction techniques under known conditions.Images were reconstructed with the Landweber simultaneous iterative reconstruction technique (SIRT) [9].

Polar Coordinate System Matrix Discretization
When populating the system matrix, the length of each ray through each pixel of the dosimeter must be calculated for every projection, since the resulting length will change as the ray is rotated in relation to the cartesian pixel grid.However, the dosimeter can be discretized into a polar pixel grid consisting of azimuthal lines and radial rings (Figure 2) [7].If the azimuthal spacing is equal to the step size between projections, the length of each ray through each pixel is constant at each projection.This allows for the calculation of ray lengths through each pixel of just a single projection thus, reducing the required size of the system matrix by a factor of the number of projections.The radial rings were asymmetrically spaced such that the area of each polar pixel was constant ensuring each pixel was equally weighted during iterative reconstruction.For display purposes, the polar pixel images are converted back to Cartesian pixel images by means of oversampling.A 5x5 grid is created for each cartesian pixel to determine the contributions of the polar pixels to each Cartesian pixel.Each intersection of the 5x5 oversampling grid that falls within a polar pixel will give a weighting factor of one over the total number of intersections (1/25).This is used as an approximate calculation of the area of the polar pixel that overlaps the cartesian pixel (Figure 2).The value of each cartesian pixel is then calculated by the weighted sum of the contributing polar pixels.Note: Ray tracing for the polar coordinate discretized system matrix is performed with the fan beam collimated such that there are no occurrences of a single detector collecting multiple raypaths, allowing for independent analysis.

Single Detector Collection of Multiple Raypaths
As a first test case, simulations were performed with the fan beam collimated to allow rays on the outer edges of the beam span to be severely refracted and considered as crossover rays.Rays from the central region of the fan beam that follow a more direct raypath are thus considered the primary rays (Figure 1).To account for detectors receiving light from both primary and crossover raypaths, the multiple paths were treated as a single ray during the iterative reconstruction.The length through each pixel within the system matrix was then weighed based on the relative intensities of the rays to account for their relative contribution to the cumulative light intensity measured by the detectors.Synthetic sinogram data that included these crossover rays was generated for a simple ring phantom.Images were iteratively reconstructed with system matrices that accounted for these crossover rays and compared to images that were reconstructed when only primary rays were taken into account.

Polar Coordinate System Matrix Discretization
Tracing 640 rays over 360 projections on a 512x512 Cartesian coordinate grid resulted in a system matrix size of 2.16 GB when stored in the sparse matrix file format.To produce a system matrix utilizing a polar coordinate pixel grid system with the same number of pixel elements as its 512x512 Cartesian coordinate counterpart, 728 radial rings and 360 azimuthal lines were used.This yielded the expected sparse system matrix size reduction of 360 times (6.0 MB).The asymmetrical spacing of the radial rings caused the central region of the image (where radial ring spacing is the largest) to have a 50% loss in spatial resolution (0.336 mm -1 at 10% MTF) (Figure 3).Although the number of azimuthal lines is fixed to the number of projections, the number of radial rings can be increased to reduce the size of the inner most pixels and thus reduce the loss of spatial resolution in the central region of the image.This increases the number of polar pixel elements and comes at a cost of increased system matrix size.Therefore, a balance between spatial resolution of the central region of the image and system matrix size must be found.Increasing the number of radial rings to 2048 provided an increase in spatial resolution to the central region of the image (0.546 mm -1 at 10% MTF) reducing the loss to only 23% and produced a sparse system matrix size reduction of 144 times (14.6 MB) (Figure 3).

Single Detector Collection of Multiple Raypaths
The occurrence of crossover rays within the optical CT data produce artifacts in the resulting iteratively reconstructed images.Specifically, when a simple ring phantom's synthetic sinogram is reconstructed a smaller radius, less intense ring artifact is produced.By merging the crossover rays with their primary ray counterpart within the system matrix, it was shown to produce a 26% reduction in the intensity of the induced ring artifact (Figure 4).Although the ring artifact persists, this first case test problem models crossover rays that are much higher in intensity than we would expect to see from reflected rays in our optical CT system.This first case test problem should therefore be looked at as an exaggerated case of the reflected rays problem.The application of merging multiple rays within the system matrix to reflected rays within the optical CT system is still under development.However, the reduction of artifacts under the exaggerated conditions described makes that application promising.

Conclusions
A system matrix size reduction of 360 times (6.0 MB) was achieved by discretizing the image pixel elements in polar coordinates, allowing system matrix generation for only a single projection.However, this reduction in system matrix size came at the cost of decreased spatial resolution in the central region of the image.A more balanced approach found a system matrix size reduction of 144 times (14.6 MB) with only minimal reduction in spatial resolution.
Merging multiple rays within the system matrix was shown to produce a slight reduction in the intensity of image artifacts.If a more substantial correction for these crossover rays can be achieved, the correction can be applied to all sources of multiple ray paths including those from reflected rays.
These novel iterative image reconstruction methods make strides toward optical CT radiochromic gel dosimetry being a practical and accurate tool for radiation therapy treatment plan verification.

Figure 1 .
Figure1.Ray tracing simulator produced visual representation of primary (green) and crossover (red) raypaths through the prototype optical CT system of[8].

Figure 2 .
Figure 2. Polar coordinate discretization with asymmetrically spaced radial rings.Polar to Cartesian conversion is achieved through oversampling which approximates the area of polar pixel that overlaps a given Cartesian pixel.

Figure 3 .
Figure 3. Reconstructed Shepp-Logan phantom synthetic data with system matrix discretized in: A) Cartesian coordinates, B) polar coordinates with 728 radial rings and C) polar coordinates with 2048 radial rings.D) The known Shepp-Logan phantom dose distribution.

Figure 4 .
Figure 4. Reconstructed ring phantom synthetic data utilizing: A) a system matrix that only accounts for primary rays, B) a system matrix that merges crossover rays.C) Image profile comparison of the ring phantom images taken at pixel row 256.The highlighted regions show the slight reduction in intensity of the ring artifact.