Proton Computed Tomography Image Reconstruction Based on the Richardson-Lucy Algorithm

This paper proposes and validates a novel Richardson-Lucy-based iterative algorithm for proton computed tomography image reconstruction, demonstrating through Monte Carlo simulations that it achieves high spatial resolution and low relative stopping power uncertainty as a promising proof-of-concept for improving proton therapy accuracy.

The Gist

The Big Picture: Seeing Inside with "Super-Beams"

Imagine you are a doctor trying to aim a laser beam at a tumor inside a patient's body. To do this safely, you need a perfect map of the body's internal landscape. You need to know exactly how "thick" or "dense" different tissues are so the laser stops exactly at the tumor and doesn't burn healthy organs behind it.

Currently, doctors use X-rays to make these maps. But X-rays are like using a flashlight to look through a foggy window; they give a general idea, but they aren't perfect for the specific type of laser (proton therapy) used to kill cancer.

Proton Therapy is a high-tech cancer treatment that uses heavy particles (protons) instead of light X-rays. The problem is, to use protons effectively, you need a map made with protons. This is called Proton Computed Tomography (pCT).

This paper is about inventing a new, smarter way to build that map.

---

The Problem: The "Wobbly" Proton

Imagine you are throwing a bowling ball down a hallway filled with people (the patient's body).

• X-rays are like throwing a straight arrow. It goes straight through, and you know exactly where it ended up.
• Protons are like that bowling ball. When it hits people, it doesn't just go straight; it bounces off shoulders, elbows, and knees. It wobbles and curves.

Because the protons wobble, it's very hard to figure out exactly where they went inside the body just by looking at where they entered and where they exited. If you try to draw the map based on a straight line, your map will be blurry and wrong.

The Solution: The "Richardson-Lucy" Algorithm

The authors of this paper proposed a new mathematical recipe (an algorithm) to fix this blurry map. They named it after two astronomers who originally used it to clean up blurry pictures of stars.

Think of the algorithm like a smart photo editor or a noise-canceling headphone:

1. The Guess: The computer starts with a blurry, guess-work map of the body.
2. The Check: It simulates thousands of protons flying through that guess-map.
3. The Comparison: It compares where the simulated protons should have landed versus where the real protons actually landed.
4. The Correction: If the real protons landed in a spot the computer didn't expect, the computer says, "Ah, my map was wrong here. I need to make this part of the tissue denser (or less dense)."
5. The Loop: It repeats this process millions of times, getting slightly better with every turn, until the map is sharp and accurate.

The authors call this the Richardson-Lucy method. It's special because it's very good at handling the "wobble" (scattering) of the protons without getting confused by the noise.

The Experiment: Testing the Recipe

To see if their new recipe worked, the team didn't use real patients yet. Instead, they built virtual plastic models (phantoms) that look like human tissue.

• The CTP528: A model with tiny, sharp lines to test if the image is sharp (like checking if a photo is in focus).
• The CTP404: A model with different colored blocks to test if the image is accurate (like checking if the colors are true).

They ran their new algorithm on a super-fast computer (a Graphics Processing Unit, or GPU—the same kind of chip in a video game console) and simulated millions of protons hitting these models.

The Results: A Clearer Picture

The results were very promising:

• Sharpness: With their new method, they could see very fine details (about 5 lines per centimeter). This is sharp enough to be useful for doctors.
• Accuracy: They could measure the "density" of the tissues with an error of less than 1%. This is the "gold standard" needed for safe cancer treatment.
• The Trade-off: They tested three types of "cameras" (detectors) to catch the protons.
  • The Ideal Camera: Perfect, no errors. (Result: Perfect image).
  • The Pixel Camera: A high-tech, expensive sensor. (Result: Very good image).
  • The Strip Camera: A cheaper, simpler sensor. (Result: Good image, but a bit blurrier).

Even with the cheaper, simpler sensors, the new algorithm managed to produce a high-quality image. This is huge because cheaper sensors mean the machine costs less, making the technology available to more hospitals.

Why This Matters

Think of this paper as the blueprint for a new kind of GPS for cancer doctors.

• Before: Doctors had to guess the terrain, leading to potential "traffic jams" (overdosing healthy tissue) or "getting lost" (missing the tumor).
• Now: With this new algorithm, they can build a 3D map that accounts for the "wobbly" nature of the protons.

The authors admit this is just the beginning (a "proof of concept"). They need to make it faster and test it in 3D (not just 2D slices). But they have shown that this mathematical "noise-canceling" trick works. It offers a path to making proton therapy more precise, safer, and accessible to more people fighting cancer.

In short: They found a new way to clean up the blurry, wobbly pictures taken by proton beams, turning a fuzzy guess into a sharp, reliable map for saving lives.

Technical Summary

1. Problem Statement

Proton therapy is a highly effective cancer treatment, but its precision relies heavily on accurate Relative Stopping Power (RSP) maps for dose planning. Current RSP estimation methods often suffer from uncertainties that can lead to under-dosing tumors or overdosing healthy tissue. Proton Computed Tomography (pCT) is proposed as a solution to measure RSP directly using the same proton beam used for treatment, thereby reducing uncertainties.

However, pCT reconstruction faces significant challenges:

• Physics Complexity: Unlike X-ray CT, protons undergo Multiple Coulomb Scattering (MCS), causing them to deviate from straight lines. This renders classical analytical reconstruction methods (like Filtered Backprojection) suboptimal.
• Computational Cost: Iterative methods (like Algebraic Reconstruction Technique - ART) that model curved proton paths are computationally intensive, often requiring hours on CPUs.
• Detector Limitations: Moving toward single-sided scanner designs (to reduce cost and complexity) introduces greater uncertainty in proton path estimation compared to double-sided designs, necessitating more robust reconstruction algorithms to maintain image quality.

2. Methodology

The authors propose a novel iterative image reconstruction algorithm based on the Richardson–Lucy (RL) method, adapted specifically for proton CT.

A. The Richardson–Lucy Adaptation

• Origin: The RL algorithm is traditionally used in astronomy and emission tomography for Poisson-distributed data.
• Adaptation: The authors reformulate the RL update rule for proton transmission imaging. Instead of a point-spread function (PSF), they utilize the system matrix ($A$), which represents the interaction between proton tracks and voxels.
• Update Rule: The algorithm iteratively updates the voxel RSP values ($x_j$) using the ratio of measured Water-Equivalent Path Length (WEPL) to predicted WEPL, normalized by the system matrix.
  $$x_j^{(k+1)} = x_j^{(k)} \frac{1}{\sum_i A_{ij}} \sum_i A_{ij} \frac{y_i}{h_i^{(k)}}$$
  Where $y_i$ is the measured WEPL and $h_i^{(k)}$ is the predicted WEPL at iteration $k$.
• Advantages: The method naturally enforces non-negativity, handles Poisson noise statistics, and is highly parallelizable, making it ideal for GPU acceleration.

B. Proton-Phantom Interaction Model

To address the lack of upstream detectors in single-sided designs, the authors implemented a probabilistic path estimation:

• Most Likely Path (MLP): Calculated using beam parameters and downstream tracker data.
• Probability Density: A Gaussian distribution is applied around the MLP to account for scattering within the phantom and measurement uncertainties.
• Spline Approximation: A third-order spline is used to approximate the curved trajectory inside the phantom for computational efficiency.

C. Simulation Setup

• Phantoms: Two standard phantoms were simulated using Geant4/GATE:
  • CTP528: For evaluating spatial resolution (line pairs per cm).
  • CTP404: For evaluating RSP accuracy (density resolution).
• Detector Models: Three configurations were tested:
  1. Ideal: No measurement errors.
  2. Silicon Pixel: Modeled after the Bergen pCT Collaboration design (high resolution).
  3. Silicon Strip: Modeled after the LLU/UCSC design (lower resolution).
• Data Volume: Approximately 43 million primary protons were simulated per slice (approx. 0.05 mGy dose).
• Hardware: Computations were performed on Nvidia 1080 Ti GPUs, processing data in chunks to manage memory.

3. Key Contributions

1. First Application of RL in pCT: This is the first study to apply the Richardson–Lucy algorithm to proton CT image reconstruction, demonstrating its viability as an alternative to ART-based methods.
2. Single-Sided Feasibility: The study proves that high-quality reconstruction is achievable with single-sided detector designs (which are cheaper and simpler) by using a probabilistic path model and the RL algorithm.
3. GPU Optimization: The algorithm is designed for GPU parallelization, utilizing sparse matrix-vector operations and ensemble averaging of partial results to improve signal-to-noise ratio while maintaining speed.
4. Comprehensive Benchmarking: The method was rigorously tested against standard phantoms (CTP528/404) under varying detector conditions, providing a baseline for future clinical development.

4. Results

The study achieved promising results after processing ~43 million protons:

Metric	Ideal Setup	Silicon Pixel	Silicon Strip
Spatial Resolution	4.88 lp/cm	3.11 lp/cm	2.33 lp/cm
Average RSP Error	-0.66%	-2.25%	-2.64%

• Spatial Resolution: The ideal setup achieved 4.88 line pairs/cm, reaching 97% of the limiting Nyquist frequency for a 1mm pixel grid. Even the realistic silicon strip setup (2.33 lp/cm) met minimum clinical requirements.
• RSP Accuracy: The ideal setup achieved an average RSP error of -0.66%, which is within the clinically required <1% uncertainty threshold. The realistic setups showed slightly higher errors (-2.25% to -2.64%) but remained competitive.
• Convergence: Visual satisfaction was achieved with only 10% of the data (4 million protons) in roughly 20 minutes, suggesting potential for faster clinical workflows.

5. Significance and Future Outlook

• Clinical Relevance: The results demonstrate that the Richardson–Lucy approach can produce RSP maps with sufficient accuracy and resolution for proton therapy treatment planning, potentially enabling the use of more cost-effective single-sided scanners.
• Performance vs. Cost: By achieving high resolution with single-sided designs, this method could significantly lower the barrier to entry for pCT systems in clinical centers.
• Future Work: The authors plan to extend the algorithm to 3D reconstruction, incorporate more detailed probability maps for proton paths, and further optimize processing times using advanced machine learning techniques and stronger GPU hardware.

In conclusion, this paper presents a proof-of-concept that the Richardson–Lucy algorithm, combined with probabilistic path modeling, is a robust, accurate, and computationally efficient candidate for next-generation proton CT image reconstruction.