# Developments in Mathematical Algorithms and Computational Tools for Proton
CT and Particle Therapy Treatment Planning

Yair Censor, Keith E. Schubert, and Reinhard W. Schulte, .
Preprint accepted for publication in IEEE Transactions on Radiation and Plasma
Medical Sciences.
August 16, 2021 Y. Censor is with the Department of Mathematics, The
University of Haifa Haifa, IsraelK. E. Schubert is with the Department of
Electrical and Computer Engineering, Baylor University, Waco, TX, USAR. W.
Schulte is with the Department of Basic Sciences, Loma Linda University Loma
Linda, CA, USA

###### Abstract

We summarize recent results and ongoing activities in mathematical algorithms
and computer science methods related to proton computed tomography (pCT) and
intensity-modulated particle therapy (IMPT) treatment planning. Proton therapy
necessitates a high level of delivery accuracy to exploit the selective
targeting imparted by the Bragg peak. For this purpose, pCT utilizes the
proton beam itself to create images. The technique works by sending a low-
intensity beam of protons through the patient and measuring the position,
direction, and energy loss of each exiting proton. The pCT technique allows
reconstruction of the volumetric distribution of the relative stopping power
(RSP) of the patient tissues for use in treatment planning and pre-treatment
range verification. We have investigated new ways to make the reconstruction
both efficient and accurate. Better accuracy of RSP also enables more robust
inverse approaches to IMPT. For IMPT, we developed a framework for performing
intensity-modulation of the proton pencil beams. We expect that these
developments will lead to additional project work in the years to come, which
requires a regular exchange between experts in the fields of mathematics,
computer science, and medical physics. We have initiated such an exchange by
organizing annual workshops on pCT and IMPT algorithm and technology
developments. This report is, admittedly, tilted toward our interdisciplinary
work and methods. We offer a comprehensive overview of results, problems, and
challenges in pCT and IMPT with the aim of making other scientists wanting to
tackle such issues and to strengthen their interdisciplinary collaboration by
bringing together cutting-edge know-how from medicine, computer science,
physics, and mathematics to bear on medical physics problems at hand.

###### Index Terms:

proton therapy, proton computed tomography, intensity-modulated therapy, blob
basis functions, superiorization, data partitioning, Monte Carlo simulation,
digital phantoms, motion-adapted reconstruction

##  I Introduction

Proton therapy is becoming increasingly common for cancer radiation therapy.
Protons afford tissue-sparing advantages that should be carefully tested
against the best available photon therapy techniques, i.e., intensity-
modulated radiation therapy (IMRT). However, this requires further development
of image-guidance and inverse planning techniques for proton and ion therapy.
This report summarizes some developments of advanced algorithmic and
computational methods for proton computed tomography (pCT) and for intensity-
modulated particle therapy (IMPT) that we developed by applying expertise in
optimization algorithms, computer science, and medical physics, related to
proton and ion therapy. We give references to published work when applicable
and include details about work-in-progress. The purpose of over-viewing our
work in this rapidly evolving field is to inspire other researchers to
contribute to it in an interdisciplinary fashion.

Both IMPT and pCT are “inverse problems”; they are processes of calculating
causal factors from a set of observations they produced. For example, in IMPT
or IMRT the fluence of individual pencil beams or beamlets will produce a dose
distribution and the clinical prescriptions play the role of “observations” in
this framework from which the proper fluences are calculated. Likewise, in
pCT, relative stopping power (RSP) is the property of tissues that results in
energy loss of protons and ions. Here, measured energy loss of individual
protons traversing the object are the “observations” and the RSP is calculated
from these observations, forming the pCT image. Starting from this premise, we
discuss the following topics mostly from the mathematical and computational
points of view:

  * •

The development of a feasibility-seeking approach and algorithms for treatment
planning of IMPT.

  * •

The superiorization methodology (SM) combined with the diagonally relaxed
orthogonal projection (DROP) feasibility-seeking algorithm.

  * •

The application of the SM to total variation superiorization (TVS) for pCT
reconstruction.

  * •

The use of blob basis functions and the reconstruction from pCT data acquired
in the presence of organ motion.

  * •

The derivative-free framework for the SM with component-wise perturbations.

  * •

The selection of appropriate digital anthropomorphic phantoms for exploring
pCT and IMPT.

  * •

The selection of simulation tools to generate pCT data and, as a forward
calculation method, to calculate dose contributions from proton pencil beams.

  * •

The development of tools to compare pCT and IMPT results with those of
concurrent methods.

  * •

The development of code optimized for general purpose graphics processing
units (GPGPU) that combines algorithmic advances, memory, and computational
tuning, required to achieve speedup.

  * •

The efforts to acquire pCT data with phantoms using the preclinical pCT
scanner at the Northwestern Medicine Chicago Proton Center (NMCPC).

##  II Model Formulations and Mathematical Algorithms for pCT and IMPT Inverse
Problems

In this section, we give an overview of some of the methods used in the
investigations of the topics listed above. If available, we point to
additional details about the methods in the published literature. From the
mathematical modeling and algorithmic points of view, we take recent
algorithmic developments and “translate” them to applications in pCT and IMPT
inverse problems. We demonstrate this approach by two specific novel
techniques that were investigated and have been published but can benefit from
further adaptation and refinement for these applications. These are the
superiorization methodology for pCT image reconstruction and the split inverse
problem paradigm for IMPT.

The core of the fully-discretized inverse problems of IMPT and pCT problems
discussed here consists of the constraints, which are dictated by the
underlying physical model. In the feasibility approach, we look at convex
feasibility problems (CFPs) of the form: find a vector $x^{\ast}\in
C:=\cap_{i=1}^{I}C_{i},$ where the sets $C_{i}\subseteq R^{J}$ are nonempty
closed convex subsets of the Euclidean space $R^{J}$, see, e.g., [1, 2, 3] or
[4, Chapter 5] for results and references on this broad topic.

Both IMPT inverse planning and pCT, as well as other inverse problems where
the underlying system has very large values of $I$ and $J$ and is often very
sparse, fall under this category. Under these circumstances, projection
methods are beneficial [1, 5, 6, 7, 8, 9]. These are iterative algorithms that
use projections onto sets, e.g., hyperplanes or half-spaces in the linear
cases, while relying on the general principle that when a family of, usually
closed and convex, sets is present, projections onto the given individual sets
are easier to perform than projections onto other sets (intersections, image
sets under some transformation, etc.) that are derived from the given
individual sets.111The fully-discretized formulation of the inverse problem of
IMRT treatment planning was not commonly accepted in the early days as it is
today. In 1982, Brahme, Roos, and Lax tried to solve the inverse problem of
radiation therapy treatment planning in its continuous (not fully-discretized)
formulation via integral inversion [10]. However, to be able to generate the
required inverse transform, they had to make unrealistic assumptions on the
model that rendered their analysis impractical. Altschuler and Censor proposed
in 1984 [11] a fully-discretized IMRT model. Here, “fully-discretized” means
that not only the irradiated volume is discretized into voxels but that also
the external radiation field is discretized into “rays” (“beamlets” or “pencil
beams” in today’s language). The initial conference report was followed by a
sequence of papers that established this approach [12, 13, 14]. To the best of
our knowledge, these and the independent 1990 paper of Bortfeld et al. [15]
were the first publications that suggested the fully-discretized approach to
the inverse problem of IMRT, see also [16] and the introduction of [17].

###  II-A The Superiorization Methodology

The superiorization methodology (SM) can be applied to the data of constrained
minimization (CM) problems of the form: $\mathrm{minimize}\left\\{\phi(x)\mid
x\in C\right\\},$ where $\phi:R^{J}\rightarrow R$ is a target function, and
$C\subseteq\Theta\subseteq{R^{J}}$ is a given feasible set defined by
constraints, see, e.g., [18]. It aims at finding a feasible point that is
superior (with respect to the target function value) to one returned by an
algorithm that is only feasibility-seeking. In doing so, SM is situated
between feasibility-seeking and full-fledged CM.

There are two main reasons why superiorization is beneficial: (i) for a
problem for which an exact CM algorithm has not yet been discovered, but there
are iterative feasibility-seeking methods that provide constraints-compatible
solutions, which can be turned by the superiorization methodology into methods
that will be practically useful for the target function reduction effort; and
(ii) when existing exact optimization algorithms are either very time
consuming or require too much computer space for large problems to be
processed by run-of-the-mill computers. On the other hand, space- and time-
efficient algorithms exist for constraints-compatibility-seeking, and these
can be turned into efficient algorithms for superiorization. Examples of such
situations are given in [19, 20].

There is no general answer yet to the question under which circumstances one
should resort to superiorization as the method of choice. Even when tractable
constrained optimization algorithms are available, the SM often yields better
or comparable results. This has been demonstrated in [19] where a comparison
between the projected subgradient method and SM showed that the SM performed
better when applying it to an image reconstruction from projections problem.
Many works cited on [21] attest to the practical success of the SM in a
variety of situations.

In the SM, one associates with the feasible set $C$ a proximity function
${Prox}_{C}:\Theta\rightarrow R_{+}$, which is an indicator of how
incompatible a vector $x\in\Theta$ is with the constraints. For any given
$\varepsilon>0$, a point $x\in\Theta$ for which ${Prox}_{C}(x)\leq\varepsilon$
is called an $\varepsilon$-compatible point for $C$.

The $\varepsilon$-output of a sequence, defined in [22], can be explained as
follows. For a nonnegative $\varepsilon$ and a sequence
$R:=\left(x^{k}\right)_{k=0}^{\infty}$ of points in $\Theta$, the
$\varepsilon$-output of the sequence $R$ is a point $x\in\Theta$ that has the
following properties: ${Prox}_{C}(x)\leq\varepsilon,$__ and there is a
nonnegative integer $K$ such that $x^{K}=x$ and, for all nonnegative integers
$k<K$, we have ${Prox}_{C}(x^{k})>\varepsilon$__. If there is such an $x$,
then it is unique. If there is no such $x$, then we say that the
$\varepsilon$-output of the sequence is undefined and otherwise, that it is
defined.

In order to “superiorize” a feasibility-seeking algorithm, commonly called the
“Basic Algorithm”, represented by the iteration
$x^{k+1}=\boldsymbol{A}_{C}(x^{k}),$ we need it to have strong perturbation
resilience in the sense that for every $\varepsilon>0,$ for which an
$\varepsilon$-output is defined for a sequence generated by the Basic
Algorithm for every $x^{0}\in\Theta$, we have also that the
$\varepsilon^{\prime}$-output is defined for every
$\varepsilon^{\prime}>\varepsilon$ and for every sequence
$\left\\{y^{k}\right\\}_{k=0}^{\infty}$ generated by

$y^{0}=x^{0},\>\>y^{k+1}=\boldsymbol{A}_{C}\left(y^{k}+\beta_{k}v^{k}\right),\text{
for all }k\geq 0,$ (1)

where the vector sequence $\left\\{v^{k}\right\\}_{k=0}^{\infty}$ is bounded
and the scalars $\left\\{\beta_{k}\right\\}_{k=0}^{\infty}$ are such that
$\beta_{k}\geq 0$, for all $k\geq 0,$ and
$\sum_{k=0}^{\infty}\beta_{k}<+\infty$. See, e.g., [19, 22] for additional
details.

Sufficient conditions for strong perturbation resilience of a Basic Algorithm
were derived in [18, Theorem 1] and [19]. Along with the constraints
$C\subseteq R^{J}$, we look at a target function $\phi:\Theta\subseteq
R^{J}\rightarrow R$, with the convention that a point in $R^{J}$ for which the
value of $\phi$ is smaller is considered superior to a point in $R^{J}$ for
which the value of $\phi$ is larger. The essential idea of the superiorization
methodology is to make use of the perturbations of (1) to transform a strongly
perturbation resilient algorithm that seeks feasibility into an algorithm with
an output that is equally good for constraints-compatibility but is also
superior (not necessarily optimal) according to the target function $\phi$.

The SM accomplishes its goal by producing from the Basic Algorithm another
algorithm, called its superiorized version, that implements, at every
iteration, a perturbation that reduces the target function value locally,
i.e., $\phi\left(y^{k}+\beta_{k}v^{k}\right)\leq\phi\left(y^{k}\right)$. The
Superiorized Version of the Basic Algorithm assumes that we have available a
summable sequence $\left\\{\eta_{\ell}\right\\}_{\ell=0}^{\infty}$ of positive
real numbers (for example, $\eta_{\ell}=a^{\ell}$, where $0<a<1$) and it
generates, simultaneously with the sequence
$\left\\{y^{k}\right\\}_{k=0}^{\infty}$ in $\Theta$, sequences
$\left\\{v^{k}\right\\}_{k=0}^{\infty}$ and
$\left\\{\beta_{k}\right\\}_{k=0}^{\infty}$. The latter is generated as a
subsequence of $\left\\{\eta_{\ell}\right\\}_{\ell=0}^{\infty}$, resulting in
a nonnegative summable sequence $\left\\{\beta_{k}\right\\}_{k=0}^{\infty}$.
The algorithm further depends on a specified initial point $y^{0}\in\Theta$
and a positive integer $N$. It makes use of a logical variable called loop _._
Such a superiorized algorithm is presented here by its pseudo-code in
Algorithm 1.

Algorithm 1: The Superiorized Version of the Basic Algorithm

set $k=0$

set $y^{k}=y^{0}$

set $\ell=-1$

repeat

set $n=0$

set $y^{k,n}=y^{k}$

while $n$$<$$N$

set $v^{k,n}$ to be a non-ascending vector for $\phi$ at $y^{k,n}$

set _loop=true_

while _loop_

set $\ell=\ell+1$

set $\beta_{k,n}=\eta_{\ell}$

set $z=y^{k,n}+\beta_{k,n}v^{k,n}$

if $\phi\left(z\right)$$\leq$ $\phi\left(y^{k}\right)$ then

set $n$$=$$n+1$

set $y^{k,n}$$=$$z$

set _loop = false_

set $y^{k+1}$$=$$\boldsymbol{A}_{C}\left(y^{k,N}\right)$

set $k=k+1$

In general, the Superiorized Version of the Basic Algorithm outputs solutions
that are essentially as constraints-compatible as those produced by the
original (not superiorized) Basic Algorithm. However, due to the repeated
steering of the process toward reducing the value of the target function
$\phi$, we can expect that the output of the Superiorized Version will be
superior (from the point of view of $\phi$) to the output of the original
algorithm.

In recent years, the algorithmic structure of the superiorization method has
undergone some evolution in ways that offer benefits in pCT. The details of
this evolution can be found in the Appendix of [23], titled “The algorithmic
evolution of superiorization”.

A comprehensive overview of the state of the art and current research on
superiorization appears in our continuously updated bibliography Internet page
that currently contains 138 items [21]. Research works in this bibliography
include a variety of reports ranging from new applications to new mathematical
results on the foundations of superiorization. A special issue entitled:
“Superiorization: Theory and Applications” of the journal Inverse Problems
[24] contains several interesting papers on the theory and practice of SM. A
special issue entitled: “Superiorization versus Constrained Optimization:
Analysis and Applications” of the Journal of Applied and Numerical
Optimization (JANO) appeared in 2020 [25].

A word About the history: The superiorization method was born when the terms
and notions “superiorization” and “perturbation resilience”, in the present
context, first appeared in the 2009 paper of Davidi, Herman and Censor [20]
which followed its 2007 forerunner by Butnariu et al. [26]. The ideas have
some of their roots in the 2006 and 2008 papers of Butnariu et al. [27] and
[28]. All these culminated in Ran Davidi’s 2010 PhD dissertation [29] and the
many papers since then cited in [21].

###  II-B Derivative-free Superiorization (DFS) with Component-wise
Perturbations

Superiorization reduces, but not necessarily minimizes, the value of a target
function while seeking constraints-compatibility. When the perturbation steps
are computationally efficient, the superior result is obtained with
essentially the same computational cost as that of the original feasibility-
seeking algorithm.

In the literature on superiorization, the perturbations that interlace target
function reduction steps into the basic algorithm have, up to now, been done
mostly by using negative gradients (or subgradients) directions and, thus,
require some form of differentiability of the target function. In [30], we
introduced component-wise perturbations that allow local non-ascent of the
target function. These enable the SM to be applied with target functions that
are not differentiable, similarly to such situations in optimization theory,
where coordinate descent approaches are used.

The ramification of DFS for practical applications of the SM to IMPT treatment
planning are meaningful. For example, normal tissue complication probability
(NTCP), which is a predictor of radiobiological effects for organs at risk,
should be used as an objective function in the mathematical problem modeling
and the planning algorithm. Using NTCP or similar functions is hampered
because these functions are, in general, empirical functions whose derivatives
cannot be calculated, see, e.g., [31]. In a recent paper [32] the authors say
that “…practical tools to handle the variable biological efficiency in Proton
Therapy are urgently demanded…” highlighting this need for the near future,
but also in a longer perspective, of the proton therapy community.

By considering component-wise perturbations, we generalized previous
superiorization schemes to enable the use of a more extensive selection of
methods for step-wise reduction of the target function. As a first step in
validating component-wise perturbations, we presented in [30] a new
superiorization scheme for reducing total variation (TV) during image
reconstruction from projections. More recently, this work has been continued
in [33]. In that paper, DFS is put in context with the large field of
derivative-free optimization (DFO) with many relevant references and a tool,
called a “proximity-target curve”, for deciding which of two iterative methods
is “better” for solving a particular problem is developed. It is worthwhile to
mention that other methodologies such as averaged stochastic gradient descent
and automatic differentiation, that are widely used in optimization theory, in
particular in the field of machine learning, also cope with lack of
differentiability and employ “coordinate-wise” searches. There is certainly a
potential to learn from the advances in these fields toward improving DFS.

###  II-C The Fully-Discretized Problem Formulation of pCT Image
Reconstruction

The water equivalent path length (WEPL) of a proton through an object, i.e.,
the length of the path the proton travels through water that leads to the same
mean energy loss as in the object, can be expressed as a line integral of
relative (to water) stopping power (RSP) along the path [34]. In practice, the
object is described by a set of basis functions, and the integrals are
expressed as discrete sums. This leads to a linear system $Ax=b,$ where the
$ij$th entry of the matrix $A$ is the intersection length of the path of the
$i$th proton with the $j$th voxel (basis function) and the $i$th component of
the measurement vector $b$ is the $i$th proton’s measured WEPL. The
mathematical formulation of the fully-discretized pCT reconstruction problem
is: Given $A$ and $b$, estimate $x$.

The system matrix $A$ is generated by calculating the path length through each
basis function (normally a voxel). Different from the straight-line assumption
in x-ray CT, in pCT, discrete steps along most likely path (MLP)
approximations developed in [35], created by multiple small-angle Coulomb-
scattering of protons in the object, are used. Thus, proton paths through the
reconstructed object are approximated at a finite set of points along each
MLP. Employing a linear path approximation through individual basis functions
simplifies the calculation. It is justified by the minor deviation of a proton
path from a straight line on the scale of the voxel dimension ($1$-$2$ mm).
However, the need to identify intersection lengths of protons through
individual voxels creates high computational demand with the need to identify
intersection lengths of millions of protons through millions of voxels. We
addressed this problem by using GPGPUs (general-purpose graphics processing
units) that perform calculations that would otherwise typically be conducted
by the CPU (central processing unit). Such GPGPUs were used for the
implementations of pCT reconstruction codes (see Subsection III-A below).

###  II-D A New Algorithm for TV-Superiorized (TVS) pCT Image Reconstruction

The efficacy of the SM for image reconstruction in pCT has been shown in
previous work [36]. In that work, we superiorized the total variation (TV) as
the target function to improve pCT image quality. The usefulness of TVS was
demonstrated for pCT image reconstruction with different superiorized versions
of the block-iterative diagonally relaxed orthogonal projections (DROP)
algorithm [37]. Two TVS schemes added-on to DROP were investigated; the first
carried out the superiorization steps once per cycle and the second once per
block.

A new version of the TVS algorithm, referred to as NTVS (New TVS), was
investigated in [38]. Compared to the original TVS algorithm published by
Penfold et al. in [36], the NTVS includes several structural changes and new
aspects previously not investigated. It combines properties that were
scattered among previous works on TVS in x-ray CT but were never combined in a
single algorithm, neither for x-ray CT nor for pCT. These properties are: (1)
exclusion of the TV reduction verification step; (2) usage of powers of a
perturbation kernel $\alpha$ to control the step-sizes $\beta_{k}$ in the TV
perturbation steps; (3) incorporation of the user-chosen integer $N$ that
specifies the number of TV perturbation steps between consecutive feasibility-
seeking iterations; and (4) incorporation of a new formula for calculating the
power of $\alpha$ that is used to calculate the step-size of the perturbation
steps. The notation and other algorithmic details of the NTVS algorithm can be
found in [38].

###  II-E Adding Robustness to pCT Image Reconstruction

Imaging systems commonly contain some uncertainty in modeling and measurement.
For pCT, there is uncertainty in the estimated MLP and the measurement of the
WEPL of individual protons. In standard algorithms, this uncertainty, often
measured by the residuals $r=\parallel b-Ax\parallel$, is considered part of
the measurements $b$ rather than the modeled system $Ax$. In robust systems,
errors are modeled as part of both the system matrix $A$ (in pCT, intersection
length of individual proton MLPs with individual voxels) and the measurement
vector $b$ (in pCT, the WEPL). Methods such as total least squares, ridge
regression, and Tikhonov regularization are classic techniques for robust
systems, and all can be solved as a diagonal perturbation to the normal
equations of least squares, see for example [39]. The main difference relates
to the regularization parameter, which is typically a small number related to
the uncertainty parameter, and bounds on $A$ and $b$. A more recent technique,
Bounded Data Uncertainty (BDU) models the uncertainty as perturbations $E$, to
the matrix $A$, such that the actual system is $(A+E)x$. The perturbations are
bounded such that $\|E\|\leq\eta$. For the worst case perturbations [40] and
[41] and best case perturbations [42] and [43] the resulting problem is solved
by $x=(A^{t}\Phi A+\Psi)^{-1}A^{T}\Phi b$, which again can be solved by a
diagonal perturbation to the normal equation of least squares. In BDU, the
perturbation to the diagonal is dependent on the solution and must be solved
for iteratively using a “secular equation,” which is a one dimensional,
nonlinear equation involving uncertainty parameter $\eta$, and bounds on $A$
and $b$. These diagonal perturbations can be thought of as perturbing the
singular values of the underlying system of $A$ and $b$. Analogously, in [44],
we introduced the _fully-simultaneous adaptive iterative solver_ (FSAIS) for
the fully-discretized problem formulation of pCT image reconstruction and
provided a convergence analysis. To obtain a sparse robust solution, we formed
the augmented linear system

$\left[\begin{array}[]{cc}-\varPsi&A^{T}\\\
A&\varPhi^{-1}\end{array}\right]\left[\begin{array}[]{c}x\\\
r\end{array}\right]=\left[\begin{array}[]{c}0\\\ b\end{array}\right],$ (2)

where $A$, $x$, and $b$ are defined as in II-C, and the choice of the matrix
parameters $\varPsi$ and $\varPhi$ is discussed below. Solving the top
equation for $x$ and the second equation for $r$, one obtains the following
system of equations:

$r=\varPhi\left(b-Ax\right),$ (3)

$x=\varPsi^{-1}A^{T}r,$ (4)

which reduces to

$x=\varPsi^{-1}A^{T}\varPhi\left(b-Ax\right).$ (5)

This equation opens a formal way of selecting an iterative algorithm that
leads to robust RSP values. Given a current iteration vector $x^{k}$, the
FSAIS generates the next iteration vector $x^{k+1}$ by

$x^{k+1}=x^{k}+\varPsi^{-1}A^{T}\varPhi(b-Ax^{k}).$ (6)

Since many robust methods result in diagonal perturbations with values that
typically lie in a similar area, in FSAIS the two diagonal parameter matrices
$\varPsi$ and $\varPhi$ are free parameters that can be adjusted for
robustness and speed of convergence. The parameter $\Phi$ handles
normalization and weighting between the rows of $A$ and $b$. The parameter
$\Psi$ includes the relaxation paramter and adjusts weighting between the
uncertainty $\eta$, the elements of $x$, and the matrix $A$. Our choices for
these parameters are $\varPhi_{i}:=\frac{{1}}{\|A_{i}\|^{2}}$, where $A_{i}$
is the $i$th row of the matrix $A$, and
$\varPsi_{j}^{-1}(k):=(1-x_{j}^{k})\lambda(k)$, where $\lambda(k)$ is the
relaxation parameter at the $k$th iteration and $x_{j}^{k}$ is the $j$th
component of the $k$th iteration vector. In a pCT reconstruction study with
the Catphan® CTP404 sensitometry module (The Phantom Laboratory, Inc., Salem,
NY), which contains cylindrical inserts with a wide range of RSP, FSAIS was
shown to obtain results similar to the DROP algorithm of [37] for RSP values
near $1$ but obtained significantly better results than DROP for denser
materials (RSP>1) such as acrylic (1.160), Delrin® (1.359), and Teflon®
(1.790). FSAIS was also able to handle missing data, i.e., proton histories in
an angle interval being removed. This was done to simulate situations, for
example, where protons cannot penetrate individual body sections because the
proton energy would not be sufficient. Using the filtered back projection
(FBP) reconstruction along straight proton paths as the initial iterate, FSAIS
stayed within 2% of the RSP for up to 60 degrees of lost data and was within
1% for almost all materials for 30 degrees of lost data.

###  II-F Development of a Motion-Adapted Iterative Reconstruction Algorithmic
Framework for pCT

In the past and up to now, it was assumed that the object that undergoes pCT
imaging is static throughout the pCT data acquisition. However, this is not
the case when imaging, for example, the lungs because, during the image
acquisition, the acquired data belong to different motion states at different
specific time instances. It is important to note here that we register
individual protons with a time tag during the data acquisition process. During
the data acquisition, one can also measure a surrogate breathing signal that
gives information about the tidal volume of the lungs. Common concepts in
four-dimensional-CT (4D-CT) imaging with the registration of the breathing
motion are those of a “reference state” and the “reference image” which
represent the patient lungs in a specific state of filling and, typically, at
rest, respectively, for example, at the end of normal expiration when the
breathing airflow ceases for a moment. In the work-in-progress described next,
we developed and investigated an algorithmic framework that enables the
reconstruction of the reference image of a breathing patient assuming that the
tissue motion can be described by accurate deformation vector fields, see,
e.g., [45]. For ordinary (motion-free) pCT data acquisition, one uses a
detector-fixed coordinate system for tracking discrete proton coordinates for
MLP calculation and an object-fixed coordinate system for reconstructing the
RSP; both coordinate systems are represented by a regularly spaced, three-
dimensional voxel grid [46]. When motion is present, the object-fixed voxel
grid is only regular for the reference state. However, it is deformed for all
other motion states according to the deformation described by the motion
model. The general reconstruction equation $Ax=b$ is retained because motion
only influences the $A$-matrix elements according to the MLP intersection
length of the deformed grid cells. Assuming that the RSP values of the
deformed grid cells are retained (an approximation), the reconstructed image
corresponds to the reference image when the solution vector is assigned to the
rectangular voxel grid, and any other motion state when it is assigned to the
deformed grid of that motion state.

We define a “motion model,” in the context of pCT image reconstruction, to be
a system (like a blackbox) capable of describing the three-dimensional (3D)
spatial displacement of human moving tissues at any time instance relative to
their reference position. To make this statement more concrete, we consider a
3D regular voxel grid fixed in the object-fixed coordinate system whose voxel
centers are, for each voxel $j$, at the 3D position vector
$\overrightarrow{X}_{0}^{(j)}.$ For every voxel index, $j=1,2,\ldots,J,$ the
pair $(y_{0}^{(j)},\overrightarrow{X}_{0}^{(j)})$ represents the relative
stopping power (RSP) of the voxel $j$, denoted by $y_{0}^{(j)},$ and the
location $\overrightarrow{X}_{0}^{(j)}$ of the voxel center in the reference
image, see Figure II.1.

Figure II.1: Two-dimensional schematic of the reference image (left) with a
regular system of points marking the centers of an object-centered fixed voxel
grid $j=1,2,\ldots,J$. In the moving image (right), the position of the
original voxel center points have shifted as a whole (translation) and
relative to each other (deformation) to new locations due to tissue motion,
which is described by a deformation vector field $T$. The blue-shaded regions
are Voronoi cells (see text for details).

According to the modeled displacement of the object tissues due to motion, at
time $t$, the pair $(y_{0}^{(j)},\overrightarrow{X}_{0}^{(j)})$ becomes
$(y_{0}^{(j)},\overrightarrow{X}_{t}^{(j)}),$ meaning that the RSP value
$y_{0}^{(j)}$ remains unchanged, but its reference point has moved to a new
spatial location $\overrightarrow{X}_{t}^{(j)},$ not necessarily conforming
with the original fixed grid, as determined by the nature of a transformation
$T$ described by a deformation vector field that is, in general, not
repeatable nor periodic and is described by the motion model.

Our motion-adapted pCT reconstruction makes use of the concept of Voronoi
diagrams and cells, see, e.g., [47, Subsection 5.5], [48]. We use the
following definition.

######  Definition 1.

Let $P:={p_{1},p_{2},\ldots,p_{J}}$ be a set of $J$ voxel center points
(called sites) in the 3D reference image space and let
$\text{dist}(p_{j},p_{k})$ denote the Euclidean distance between voxel center
point $p_{j}$ and any other voxel center point $p_{k},k\neq j$. The “Voronoi
diagram” of $P$, is defined as the partition of the 3D space into $J$
connected 3D cells (called “Voronoi cells” or “Voronoi polygons”) such that
any point $q$ in the image space belongs to the Voronoi cell corresponding to
a site $p_{j}$ if and only if $\text{dist}(q,p_{j})<\text{dist}(q,p_{k})$ for
each $p_{k}\in P$ with $k\neq j$.

It is understandable from Figure II.1, when expanded into the third spatial
dimension, that boundaries of Voronoi cells are formed by plane sections that
are created by the intersection of bisecting half-spaces between pairs of
sites, line segments formed by the intersection of these planes, and vertices
formed by the intersection of two or more of these line segments. The
“bisecting plane” of two points $p$ and $q$ is the plane that is perpendicular
to the line segment $\overline{pq}$ and intersects it at the midpoint. This
bisecting plane divides the 3D space into two half-spaces. Then the Voronoi
cell of point $p_{k}$, denoted by $\mathcal{V}(p_{k})$, is

$\mathcal{V}(p_{k})=\bigcap_{j\neq k,j=1}^{J}h(p_{k},p_{j}),$ (7)

where $h(p_{k},p_{j})$ is the bisecting plane of the points $(p_{k},p_{j})$.

Figure II.2: Two-dimensional schematic of the pCT reconstruction problem of a
moving organ at a specific time instance $t$ when the organ is represented by
the Voronoi cell basis functions shown. The most likely path (MLP) of an
individual proton intersects the shaded cells. The summation of the
intersection lengths of the $i$th MLP with the $j$th voxel, $a_{ij}$
multiplied by the unknown RSP values $y_{0}^{j}$ of each cell is set equal to
the $i$th component of the $b$ vector of measured WELP values, and thus form a
linear system of equations $Ay_{0}=b$ that is solved iteratively for the
reference image vector $y_{0}$ by a feasibility-seeking algorithm.

In our approach to motion-adapted pCT reconstruction, which is illustrated in
Figure II.2, the Voronoi cells become the new basis functions, and the
originally formulated fully-discretized problem of pCT reconstruction remains
unchanged. The specific Voronoi cell basis functions depend on the motion
state of the tissues, which is described by a motion model. The motion model
allows calculation of the displacement vector for individual voxel points at
any given instance of time. The new positions of the voxel points are used to
compute the plane segments of the Voronoi cells. Together with an MLP
algorithm, one can then calculate the intersection length of the $i$th proton
through the $j$th Voronoi cell basis function, and the $i$th component of the
WEPL measurement vector $b$ gives the $i$th proton’s measured WEPL. This
results in a linear system of equations that is solved iteratively by a
feasibility-seeking algorithm. The unknown RSP values of each cell are assumed
not to change from their values in the reference image. Thus, the
reconstructed voxel image where each voxel is assigned the reconstructed RSP
value $y_{0}^{(j)},j=1,2,\ldots,J,$ corresponds to the reference image. One
can also calculate the image at any point in time by calculating a weighted
average of the RSPs with weights that are proportional to the area a voxel
shares with overlapping Voronoi cells.

The motion model parameters could be known at the time of pCT imaging from a
separate motion imaging study, or they could be included in the image
reconstruction problem and estimated after the imaging. By way of example
only, we mention here the lung motion model originally developed by Low and
colleagues [49, 50]. This model has recently been applied in conjunction with
an iterative projection algorithm for a motion-compensated simultaneous
algebraic reconstruction technique (MC-SART) of cone beam CT (CBCT) [51, 52].
There is, however, a fundamental conceptual difference between our approach
and the MC-SART algorithm. To “handle” the motion, the authors modified the
reconstruction algorithm (SART in their work) while we, to handle the motion,
modify the problem (the linear system of equations in our work). Our approach
is, therefore, general in the sense that once the problem is appropriately
modified to handle the motion, one can apply any iterative algorithm to the
“modified problem.” Our framework can take the output of any motion model as
input to the algorithm for computation of the basis functions at each instance
in time when proton histories are acquired. The motion model and the algorithm
are, thus, two separate entities that have to be used in tandem when applied
to simulated or real data.

Besides the motion model of Low et al., other motion models exist, see, e.g.,
[53, 54, 55].

###  II-G Blob Basis Functions for pCT Image Reconstruction

“Blob” basis functions, called “radial basis functions” in approximation
theory, see, e.g., [56], are spherically symmetric basis functions for CT
reconstruction that were suggested in the field of image reconstruction from
projections by Lewitt in the early 1990s [57, 58]. These blobs are
generalizations of a well-known class of functions used in digital signal
processing called Kaiser-Bessel window functions and they yield excellent
results in multiple imaging modalities, see, e.g., the enlightening tutorial
by Herman [59]. Matej and Lewitt [60, 61] provided a careful investigation of
how the blob basis functions should be chosen when they are used in the
context of image reconstruction from projections. Since then blobs have been
used extensively for image reconstruction in X-ray computerized tomography,
positron emission tomography, single photon emission computerized tomography,
optoacoustic tomography and electron microscopy (consult [59] for references
and details). Adhering to the term “blobs”, as commonly used in this field, we
have started to develop a methodology for using blobs instead of standard
square voxels in pCT. The use of blobs will potentially improve pCT image
quality and computational efficiency.

Current pCT makes use of conventional voxels for reconstructing and
representing reconstructed images. However, proton paths through the
reconstruction object are curved and must be approximated at a finite set of
points along each path. The computational burden to identify intersection
lengths of protons through individual voxels is quite large, especially when
the object elements are not spatially invariant. Using blobs as a replacement
of voxels in pCT can significantly reduce this computational burden because of
the spherical symmetry of the blobs. We presented initial encouraging results
in collaboration with UCLA graduate student Howard Heaton as a poster at the
2016 Joint Mathematics Meetings (JMM) [62].

Figure II.3: Profile of a voxel (a) and a blob (b).

Voxels have a uniform value inside a set domain, and blobs are spherically
symmetric and taper smoothly to zero at their border, see Figure II.3. The
formal definition for blob basis functions for CT imaging, first proposed by
Lewitt [58], is as follows:

$b_{j}(r):=\begin{cases}\frac{\big{[}1-(r/a)^{2}\big{]}^{m/2}}{I_{m}(\alpha)}I_{m}\big{(}\alpha\sqrt{1-(r/a)^{2}}\big{)},&\text{if
$r\in[0,a],$}\\\ 0,&\text{otherwise},\end{cases}$ (8)

where $b_{j}$ is the basis function of the $j$th blob, $r$ is the radial
distance from the center of the $j$th blob, $I_{m}$ denotes the modified
Bessel function of the first kind of order $m$, $a$ denotes the radius of the
blobs, and $\alpha$ is a nonnegative real number that controls the shape and
taper of the blob. It is common practice to choose $m=2$, which gives a smooth
and differentiable blob function. Following Benkarroum et al. [63], we have
chosen the blob parameters to be $a=2.453144$ mm and $\alpha=13.738507$. Those
parameter values were chosen for optimal representation of piecewise-constant
images.

To generate the system matrix $A$, the path length through each basis function
must be computed. Employing a linear path approximation through individual
basis functions (voxels or blobs) allows the use of the Radon transform
$\mathcal{R}$ to approximate each entry ${A}_{ij}$, i.e.,

${A}_{ij}=\int_{path}b_{j}(\vec{s}\hskip 5.0pt)\ d\vec{s}\ \approx\
[{\mathcal{R}}b_{j}](\ell_{i},\theta_{i}),$ (9)

where

$[\mathcal{R}b_{j}](\ell_{i},\theta_{i})=\int_{line}b_{j}\left(\sqrt{\ell_{i}^{2}+z^{2}},\
\theta_{i}+\tan^{-1}\left({z}/{\ell_{i}}\right)\right)\ dz$ (10)

is the Radon transform of the straight-line approximation of the path through
the blob. Note that with blobs, the Radon transform depends solely on
perpendicular distance from the blob center, $\ell_{i}$.

To construct the system matrix, we use the following problem formulation.

######  Problem 1.

Generation of the $A$-matrix for blob basis functions: Let
$\\{\vec{\varphi}^{\hskip 2.5pt\ell}\\}_{\ell=1}^{L}$ denote an ordered set of
points along a proton path in $\mathbb{R}^{3}$. Suppose that this path passes
through an object represented with a set $\\{b_{j}\\}$ of blob basis
functions. Using successive points $\vec{\varphi}^{\hskip 2.5pt\ell}$ along
the path, uniquely estimate each nonzero blob intersection length of the path
to generate the corresponding system matrix.

A key step in the algorithm presented in [62] is to identify blobs within a
given proximity of a point $\vec{\varphi}^{\hskip 2.5pt\ell}$ on the path.
This is accomplished by identifying a corresponding point
$\vec{\varPhi}^{\hskip 2.5pt\ell}$ located nearby on the grid of blob centers
with grid unit $\beta.$ The algorithm identifies each nonzero blob
intersection length by identifying a blob in proximity to
$\vec{\varphi}^{\hskip 2.5pt\ell}$. Then the algorithm begins to cycle through
nearby blobs within a chosen range $\kappa$. Finally, it restricts
intersection lengths to being assigned during the last step the proton takes
before passing the blob center in depth. Further details can be found in [62].

###  II-H Algorithms for IMPT: Dose-volume Constrained Split Feasibility

IMPT is a rapidly developing field wherein active pencil beam scanning became
the norm in proton and ion beam treatment centers. IMPT requires scanning of a
particle treatment field with a grid of pencil beam spots that are rapidly and
often repeatedly “visited” during the treatment to deliver a spot dose with
varying intensity. Such a scan demands fast IMPT algorithms that solve the
inverse problem with multiple dose constraints within target volumes and
organs at risk. During our research, we tackled the additional complication
that dose constraints are often modified, by the treatment planner, allowing
violation of the constraints by a specific percentage of the dose for a
specific percentage of the volume. The method we developed to handle this is
also useful for x-ray IMRT, see [64].

The linear feasibility problem (LFP) formulation, which is a special case of
the well-known convex feasibility problem (CFP), see, e.g., [1], forms a basic
model for the inverse problem in the fully-discretized approach to both IMPT
and IMRT treatment planning. In 2015, we showed that IMPT inverse planning is
possible by using a fully-discretized model and a feasibility-seeking
algorithmic approach [65]. In particular, we demonstrated on a simple 2D
example that solutions meeting the planning objectives could be found by these
feasibility-seeking iterative projection algorithms.

In the planning of IMPT or IMRT, one uses dose-volume constraints (DVCs) to
evaluate treatment plans. Mathematically, DVCs are percentage-violation
constraints (PVCs). PVCs single out individual subsets of the existing volume
constraints and allow a specified percentage of dose constraints to be
violated in each subset. Without incorporating the PVCs into the mathematical
inverse planning model and algorithm themselves, it is not possible to
guarantee that an appropriate solution will be found [64].

A tractable model and an algorithmic approach to solve the IMPT/IMRT inverse
planning problem as a feasibility problem that includes DVCs has been
developed. A rigorously defined notion of PVCs cause integers to enter the
problem, which makes it difficult to solve. To circumvent this difficulty, we
reformulated the PVC with the aid of a “sparsity-norm” that counts the number
of nonzero entries in a vector. There is a rich literature on how to handle
this norm, called also the $\ell_{0}$-norm, which is still stimulating recent
and ongoing related work. It is not in the scope of this paper to review this
field but see, e.g., [66] and its references, where the problem is formulated
and handled as a full fledged constrained optimization problem. Using the
$\ell_{0}$-norm enables to enforce the DVCs and leads to redefining the
“linear feasibility problem with PVCs” as another feasibility problem that
includes non-convex constraints for the sparsity-norm.

######  Problem 2.

[64] Linear Interval Feasibility with DVC for the inverse problem in the
fully-discretized approach to IMRT treatment planning. Find $x^{\ast}\in
R^{n}$ for which

$\displaystyle 0\leq A_{1}x\leq(1+\beta)b^{1},$ (11) $\displaystyle b^{3}\geq
A_{2}x\geq b^{2},$ (12) $\displaystyle 0\leq A_{3}x\leq b^{4},$ (13)
$\displaystyle x\geq 0,$ (14)
$\displaystyle\|(A_{1}x-b^{1})_{+}\|_{0}\leq\alpha m_{1},$ (15)

where $A_{1}\in R_{+}^{m_{1}\times n}$, $A_{2}\in R_{+}^{m_{2}\times n}$,
$A_{3}\in R_{+}^{m_{3}\times n}$ are given matrices, $b^{1}\in R_{+}^{m_{1}}$,
$b^{2},b^{3}\in R_{+}^{m_{2}}$, $b^{3}\in R_{+}^{m_{3}}$ are given vectors,
and $\beta>0$ and $\alpha\in[0,1]$ are given real numbers.

In this problem, applicable to both IMPT and IMRT, the sparsity constraint
(15) takes place in the space $R^{m_{1}},$ where the vectors of doses in the
organs at risk reside. Therefore, we cannot use plain feasibility-seeking
methods but need feasibility-seeking methods for “split feasibility problems.”
Thus, we recognized that Problem 2 is a split feasibility problem. Split
feasibility problems were introduced first in [67] and further studied in [68,
69] and many other publications. The reader may consult the brief review of
“split problems” formulations and solution methods in [70] for more details
and references. For the solution of Problem 2, which includes a non-convex
sparsity-norm induced constraint, we developed a new iterative projection
algorithm, which is a combination of the $CQ$-algorithm [71] and the automatic
relaxation method (ARM) [72]. Full details of this approach appear in [64],
which applied it to a single OAR. Following this line of development, we
expanded the model to include DVCs for several OARs in a more recent paper
[70]. Handling DVCs in radiation therapy treatment planning is a viable
research area and a variety of techniques have been applied to solve it, see,
e.g., [73] and references therein.

##  III Development of software tools and phantoms to test advanced pCT and
IMPT algorithms

The demand for high-performance computing capabilities for pCT image
reconstruction and IMPT algorithm development, is apparent. In this section,
we summarize developments facilitating fast algorithms for pCT and IMPT.

###  III-A GPGPU Data Partitioning

An essential aspect of our work as members of the pCT collaboration has been
to provide tools for fast and efficient image reconstruction using GPGPU
computations. The pCT detector hardware was developed in two phases from a
fairly small and slow detector system in Phase I to a much faster Phase II
system completed in 2015 [74]. As the reconstructed image gets larger, the
number of proton histories that enter the reconstruction is up to one hundred
times the number of voxels in the image, see [75]. Currently, objects
requiring a six minutes scan time (up to 360 million proton histories) with
the “Phase II proton CT scanner” can be reconstructed on the order of a few
minutes [76, 77, 78]. As we image larger objects with protons, not only does
the reconstruction time grow significantly, but, more importantly, the memory
requirements grow too. Therefore, we need to split the data into coherent data
blocks for efficient reconstruction. Since the Phase II pCT scanner acquires
the data during a 360-degree rotation of the object relative to the pCT
scanner with a wobbled or scanned beam spot [77, 79], we can efficiently split
the data into slices orthogonal to the rotational axis, see [78].

For a nominal proton entrance energy of 200 MeV and a few tens of MeV exit
energy, proton path histories within three standard deviations of the nominal
path in lateral direction do not pass through more than three 1 mm slices. So,
when reconstructing a slice, no more than three slices above or below the
slice need to be available in memory of a particular GPGPU to reconstruct the
image efficiently, see [80]. When reconstructing slice-by-slice, we can
allocate up to one GPGPU per slice. For proton histories passing through
multiple slices, each GPU associated with a slice holds the most recent copy
of the RSP values of adjacent slices. Updates from adjacent slices are
transferred when the thread calculating the update finishes an iterative cycle
of the algorithm.

The results presented in [78, 80] show that this method does not cause a
significant image quality reduction while permitting reconstruction of the
anthropomorphic head phantom described below in under 10 minutes on NVIDIA
P100 GPGPUs. Further speedup to reconstruction in under 1 minutes is within
reach.

###  III-B Monte Carlo Simulation Tools

Monte Carlo (MC) simulations are an essential tool in the development and
testing of algorithms for pCT image reconstruction and serve as a forward
calculation tool for the radiation dose deposited in object voxels of IMPT
plans. While testing algorithms on real-world experimental data is equally
important, MC-based data are more readily available and provide a great deal
of flexibility when testing and developing new algorithms or algorithmic
variants.

We have made extensive use of the general-purpose MC simulation tool Geant4
[81] to develop a faithful model of the Phase II pCT detector system developed
by the pCT collaboration (see Section III-A). For this project, we reproduced
the exact geometry and composition of materials of the scanner as well as the
research proton beamline at the Loma Linda University Medical Center. Details
of the initial Geant4 platform development and the methods we used to validate
its output, in connection with reconstruction algorithms that we developed,
can be found in Dr. Giacometti’s thesis publications, see [82]. We then used
pCT simulations produced by that platform in addition to experimental data to
test the performance of the new TVS algorithms described in Section II-A.

Later, we implemented the Phase II scanner in the TOPAS (TOol for PArticle
Simulation) platform, which is also based on the Geant4 toolkit but has a more
user-friendly interface [83]. We also tested the Phase II pCT scanner at the
Northwestern Medicine Chicago Proton Center (NMCPC), where it is operated on a
clinical proton beamline. The TOPAS model of the NMCPC beamline was validated
with beam profile measurements acquired at NMCPC and by acquiring and
simulating pCT images using a wobbled 200 MeV proton beam. With the TOPAS
tool, we then investigated the sources of systematic uncertainty introduced
during the steps of iterative RSP reconstruction with the DROP-TVS algorithm
[36] for CT QA phantom modules (Catphan®, The Phantom Laboratory Incorporated,
Salem, NY, USA) and a digital representation of a pediatric patient with known
ground truth RSP values in the simulation, see [84]. The important point here
is that with Monte Carlo simulations of the pCT system performance, we can
intentionally include or exclude different sources of random and systematic
uncertainties that arise from random fluctuations in the proton energy loss
and scattering, which are unavoidable, and sources of imperfections in the
beam source or the detector construction. For example, we had to understand
the origin of ring artifacts seen in the reconstruction of homogeneous
phantoms. Through systematic Monte Carlo studies, we learned that these
artifacts were related to using a WEPL calibration phantom with discrete
thickness steps. A simulation of the pCT detector replacing the step phantom
with a wedge phantom demonstrated that the ring artifacts could be
significantly reduced [84].

At the end of 2016, we shipped the Phase II pCT scanner to the Heidelberg Ion
Therapy (HIT) Center in Germany for first experiments with helium ion imaging.
Again, TOPAS was used to simulate helium CT data for predicting and evaluating
the characteristics of HeCT images reconstructed with the DROP-TVS algorithm
now in comparison with proton CT images, see [85]. Given the ease and accuracy
of results, it is thus a “natural” choice to use TOPAS for forward dose
calculations when testing new IMPT algorithms.

###  III-C Digital Phantoms

Another critical need is to identify or develop anthropomorphic digital
phantoms necessary for realistic pCT and IMPT algorithm testing. The general
strategy we have followed for most of the algorithms we have tested has been
to demonstrate the algorithmic idea with several phantoms ranging from simple
phantoms often in two dimensions to increasingly larger and more realistic
phantoms. We believe that this approach prevents the “masking” of the
principal numerical performance of new algorithms by difficulties stemming
from the complexity of realistic object data, which may require different
adjustments after the performance of the algorithms has been demonstrated on
simpler datasets. In the following paragraph, we describe the digital head
phantom that we specifically designed to support our pCT algorithmic
developments.

#### Head Phantom

We created the HIGH_RES_HEAD digital head phantom with very high resolution,
[86] for the development and evaluation of pCT image reconstruction algorithms
and IMPT algorithms, but it can be used for other medical physics purposes as
well. The phantom is now available as a computational DICOM head phantom to
the Geant4 user community in the open source Geant4 MC simulation tool package
[87]. The phantom represents a head of a 5-year-old child. It is a very close
approximation of the commercial head phantom (model HN715, CIRS, Norfolk, VA,
USA). The HIGH_RES_HEAD digital phantom includes most anatomical details of a
human head and is characterized by a high spatial resolution ($0.18\times
0.18\times 1.25~{}mm^{3}$ voxel size). We used it together with the Geant4
simulation platform of the pCT system to validate its performance by comparing
digital phantom reconstructions with reconstructions from experimental data
[88]. The phantom was also used for testing the parameter space of a novel
superiorization algorithm for pCT reconstruction [38], described in Subsection
II-D above. Figure III.1 shows a representative proton CT reconstruction from
that work.

Figure III.1: Representative pCT reconstruction of the slice of the pediatric
head phantom containing regions of interest (left); the regions of interest
are filled in white and labeled by their composition (right). (reproduced from
[38]).

##  IV The Loma Linda University algorithm workshops

We have organized, since 2015, five “Loma Linda University (LLU) Workshops on
Algorithms and Computational Techniques in Proton Imaging and Intensity-
Modulated Proton Therapy.” For details about the previous and planned Annual
Loma Linda workshops, see http://ionimaging.org/. The first three workshops
were face-to-face meetings at LLU, Loma Linda, CA, USA, during the weeks
following the annual meetings of the American Association of Physicists in
Medicine (AAPM). Participation was by invitation only, and we welcomed applied
mathematicians, computer scientists, medical physicists, radiobiologists, and
radiation oncologists. We have organized the workshops in an informal setting:
there were no workshop fees, no parallel sessions, and no strictly timed
talks. The emphasis was on scientific exchange, learning, and discussions. In
all workshop meetings, the participants learned from each other, and new
research projects and collaborations have emerged that are still on-going. For
the 4th and 5th workshops, we also allowed participation by video-conferencing
via Zoom, and the 6th workshop in July 20-22, 2020, was a Zoom-only workshop
due to the COVID-19 pandemic. The 7th workshop was held in August 2-4, 2021 at
Loma Linda via Zoom, see: http://ionimaging.org/llu2021-overview/.

##  V Discussion and Conclusion

When it comes to advanced techniques such as pCT and IMPT, there is a strong
need for input from experts in computer science and mathematics to support the
medical physicist’s insight and understanding. With this in mind, we
collaborate on developing advanced reconstruction algorithms for pCT imaging
based on mathematical expertise in iterative projection methods. With computer
science expertise, we have achieved progress and fast implementation of novel
mathematical algorithms on advanced computing hardware. There is still much
more to accomplish with increasing demands for high-speed and high-performance
medical imaging and verification techniques leading to online-adaptive
particle therapy, the inclusion of biological weighting in treatment planning
optimization or superiorization, and a combination of different image guidance
technologies such as MRI and pCT. Such interdisciplinary collaboration should
stimulate others to follow a similar path and widen the network of
collaborations in the disciplines of medical physics, mathematical algorithms,
and computer science.

Acknowledgments. We thank the many collaborators on the projects mentioned
here, in particular, Blake Schultze, Paniz Karbasi, Caesar Ordonez, Christina
Sarosiek, Howard Heaton, Jose Ramos-Mendez, Pierluigi Piersimoni, Robert
Jones, Robert Johnson, Stuart Rowland, and Tai Dou as well as Mark Pankuch and
his medical physics and accelerator support team at the Northwestern Medicine
Chicago Proton Center for help with collecting proton CT data with the Phase
II pCT scanner. The interdisciplinary efforts described here were supported by
the United States-Israel Binational Science Foundation (BSF) ([89], [90]). The
work of Y.C. was supported by the ISF-NSFC joint research program grant No.
2874/19. All authors contributed equally to the writing of this paper. All
authors read and approved the final manuscript. We appreciate the anonymous
reviewer’s report which helped us improve the paper.

## References

  * [1] H. H. Bauschke and J. M. Borwein, “On projection algorithms for solving convex feasibility problems,” _SIAM Review_ , vol. 38, pp. 367–426, 1996.
  * [2] C. L. Byrne, _Applied Iterative Methods_. A. K. Peters Ltd., Wellesley, MA, USA, 2008.
  * [3] J. W. Chinneck, _Feasibility and Infeasibility in Optimization: Algorithms and Computational Methods_. International Series in Operations Research & Management Science, Springer Science+Business Media, LLC, 2008.
  * [4] Y. Censor and S. A. Zenios, _Parallel Optimization: Theory, Algorithms, and Applications_. Oxford University Press, New York, NY, USA, 1997.
  * [5] D. Butnariu, Y. Censor, P. Gurfil, and E. Hadar, “On the behavior of subgradient projections methods for convex feasibility problems in Euclidean spaces,” _SIAM Journal on Optimization_ , vol. 19, pp. 786–807, 2008.
  * [6] C. Byrne and Y. Censor, “Proximity function minimization using multiple Bregman projections, with applications to split feasibility and Kullback–Leibler distance minimization,” _Annals of Operations Research_ , vol. 105, pp. 77–98, 2001.
  * [7] Y. Censor, A. R. De Pierro, and M. Zaknoon, “Steered sequential projections for the inconsistent convex feasibility problem,” _Nonlinear Analysis: Theory, Methods & Applications_, vol. 59, pp. 385–405, 2004.
  * [8] Y. Censor, A. Motova, and A. Segal, “Perturbed projections and subgradient projections for the multiple–sets split feasibility problem,” _Journal of Mathematical Analysis and Applications_ , vol. 327, pp. 1244–1256, 2007.
  * [9] H. H. Bauschke and V. R. Koch, “Projection methods: Swiss army knives for solving feasibility and best approximation problems with halfspaces,” _Contemporary Mathematics_ , vol. 636, pp. 1–40, 2015.
  * [10] A. Brahme, J.-E. Roos, and I. Lax, “Solution of an integral equation encountered in rotation therapy,” _Physics in Medicine and Biology_ , vol. 27, pp. 1221–1229, 1982.
  * [11] M. D. Altschuler and Y. Censor, “Feasibility solutions in radiation therapy treatment planning,” _In: Proceedings of the Eighth International Conference on the Use of Computers in Radiation Therapy_ , pp. 220–224, 1984.
  * [12] Y. Censor, M. D. Altschuler, and W. D. Powlis, “A computational solution of the inverse problem in radiation-therapy treatment planning,” _Applied Mathematics and Computation_ , vol. 25, pp. 57–87, 1988.
  * [13] ——, “On the use of Cimmino’s simultaneous projections method for computing a solution of the inverse problem in radiation therapy treatment planning,” _Inverse Problems_ , vol. 4, pp. 607–623, 1988.
  * [14] Y. Censor, W. D. Powlis, and M. D. Altschuler, “On the fully discretized model for the inverse problem of radiation therapy treatment planning,” _in: Proceedings of the Thirteenth Annual Northeast Bioengineering Conference, (K. R. Foster, Editor), Institute of Electrical and Electronics Engineers (IEEE), New York, NY, USA_ , vol. 1, pp. 211–214, 1987.
  * [15] T. Bortfeld, J. Bürkelbach, R. Boesecke, and W. Schlegel, “Methods of image reconstruction from projections applied to conformation radiotherapy,” _Physics in Medicine & Biology_, vol. 35, pp. 1423–1434, 1990.
  * [16] Y. Censor, “Mathematical aspects of radiation therapy treatment planning: Continuous inversion versus full discretization and optimization versus feasibility,” _In: Computational Radiology and Imaging: Therapy and Diagnosis, C. Börgers and F. Natterer (Editors), The IMA Volumes in Mathematics and its Applications, Vol. 110, Springer-Verlag, New York, NY, USA_ , pp. 101–112, 1999.
  * [17] P. S. Cho and R. J. Marks II, “Hardware-sensitive optimization for intensity modulated radiotherapy,” _Physics in Medicine and Biology_ , vol. 45, pp. 429–440, 2000.
  * [18] G. T. Herman, E. Garduño, R. Davidi, and Y. Censor, “Superiorization: An optimization heuristic for medical physics,” _Medical Physics_ , vol. 39, pp. 5532–5546, 2012.
  * [19] Y. Censor, R. Davidi, G. T. Herman, R. W. Schulte, and L. Tetruashvili, “Projected subgradient minimization versus superiorization,” _Journal of Optimization Theory and Applications_ , vol. 160, pp. 730–747, 2014.
  * [20] R. Davidi, G. T. Herman, and Y. Censor, “Perturbation-resilient block-iterative projection methods with application to image reconstruction from projections,” _International Transactions in Operational Research_ , vol. 16, pp. 505–524, 2009.
  * [21] Y. Censor, “Superiorization and perturbation resilience of algorithms: A bibliography compiled and continuously updated,” _Available at: http://math.haifa.ac.il/yair/bib-superiorization-censor.html#top_.
  * [22] Y. Censor, R. Davidi, and G. T. Herman, “Perturbation resilience and superiorization of iterative algorithms,” _Inverse problems_ , vol. 26, no. 65008, 2010.
  * [23] Y. Censor, “Can linear superiorization be useful for linear optimization problems?” _Inverse Problems_ , vol. 33, p. 044006, 2017.
  * [24] Y. Censor, G. T. Herman, and M. Jiang, Eds., _Superiorization: Theory and Applications_ , vol. 33, 2017.
  * [25] A. Gibali, G. T. Herman, and C. Schnörr, Eds., _Superiorization versus Constrained Optimization: Analysis and Applications_ , vol. 2, 2020\.
  * [26] D. Butnariu, R. Davidi, G. T. Herman, and I. G. Kazantsev, “Stable convergence behavior under summable perturbations of a class of projection methods for convex feasibility and optimization problems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 1, pp. 540–547, 2007.
  * [27] D. Butnariu, S. Reich, and A. J. Zaslavski, “Convergence to fixed points of inexact orbits of Bregman monotone and of nonexpansive operators in Banach spaces,” in _Fixed Point Theory and its Applications, Conference Proceedings, Guanajuato, Mexico, 2005_ , H. F. Nathansky, B. G. de Buen, K. Goebel, W. A. Kirk, and B. Sims, Eds. Yokahama Publishers, Yokahama, Japan, 2006, pp. 11–32.
  * [28] ——, “Stable convergence theorems for infinite products and powers of nonexpansive mappings,” _Numerical Functional Analysis and Optimization_ , vol. 29, pp. 304–323, 2008.
  * [29] R. Davidi, “Algorithms for Superiorization and their Applications to Image Reconstruction,” Ph.D. dissertation, Department of Computer Science, The City University of New York, NY, USA, 2010. [Online]. Available: http://gradworks.umi.com/34/26/3426727.html
  * [30] Y. Censor, H. Heaton, and R. Schulte, “Derivative-free superiorization with component-wise perturbations,” _Numerical Algorithms_ , vol. 80, pp. 1219–1240, 2019.
  * [31] H. A. Gay and A. Niemierko, “A free program for calculating EUD-based NTCP and TCP in external beam radiotherapy,” _Physica Medica_ , vol. 23, pp. 115–125, 2007.
  * [32] H. Nystrom, M. F. Jensen, and P. W. Nystrom, “Treatment planning for proton therapy: what is needed in the next 10 years?” _The British Journal of Radiology_ , vol. 93, p. 1107, 2020.
  * [33] Y. Censor, E. Garduño, E. S. Helou, and G. T. Herman, “Derivative-free superiorization: Principle and algorithm,” _Numerical Algorithms_ , vol. 88, pp. 227–248, 2021.
  * [34] R. F. Hurley, R. W. Schulte, V. A. Bashkirov, A. J. Wroe, A. Ghebremedhin, H. F.-W. Sadrozinski, V. Rykalin, G. Coutrakon, P. Koss, and B. Patyal, “Water-equivalent path length calibration of a prototype proton CT scanner,” _Medical Physics_ , vol. 39, pp. 2438–2446, 2012.
  * [35] R. W. Schulte, S. N. Penfold, J. Tafas, and K. E. Schubert, “A maximum likelihood proton path formalism for application in proton computed tomography,” _Medical Physics_ , vol. 35, pp. 4849–4856, 2008.
  * [36] S. N. Penfold, R. W. Schulte, Y. Censor, and A. B. Rosenfeld, “Total variation superiorization schemes in proton computed tomography image reconstruction,” _Medical Physics_ , vol. 37, pp. 5887–5895, 2010.
  * [37] Y. Censor, T. Elfving, G. T. Herman, and T. Nikazad, “On diagonally-relaxed orthogonal projection methods,” _SIAM Journal on Scientific Computing_ , vol. 30, pp. 473–504, 2008.
  * [38] B. Schultze, Y. Censor, P. Karbasi, K. E. Schubert, and R. W. Schulte, “An improved method of total variation superiorization applied to reconstruction in proton computed tomography,” _IEEE Transactions on Medical Imaging_ , vol. 39, pp. 294–307, 2020.
  * [39] G. H. Golub and C. F. Van Loan, _Matrix Computations_ , 3rd ed. The Johns Hopkins University Press, Baltimore, MD, USA, 1996.
  * [40] S. Chandrasekaran, G. H. Golub, M. Gu, and A. H. Sayed, “Parameter estimation in the presence of bounded data uncertainties,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 19, pp. 235–252, 1998.
  * [41] G. A. Watson, “Data fitting problems with bounded uncertainties in the data,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 22, pp. 1274–1293, 2001.
  * [42] S. Chandrasekaran, G. H. Golub, M. Gu, and A. H. Sayed, “An efficient algorithm for a bounded errors-in-variables model,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 20, pp. 839–859, 1999.
  * [43] S. Chandrasekaran, M. Gu, A. H. Sayed, and K. E. Schubert, “The degenerate bounded errors-in-variables model,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 23, pp. 138–166, 2001.
  * [44] P. Karbasi, K. E. Schubert, B. Schultze, V. Bashkirov, R. P. Johnson, and R. W. Schulte, “Robust iterative methods: Convergence and applications to proton computed tomography,” _AIP Conference Proceedings_ , vol. 2160, 050008, 2019\. [Online]. Available: https://doi.org/10.1063/1.5127700.
  * [45] D. Thomas, J. Lamb, B. White, S. Jani, S. Gaudio, P. Lee, D. Ruan, M. McNitt-Gray, and D. Low, “A novel fast helical 4D-CT acquisition technique to generate low-noise sorting artifact-free images at user-selected breathing phases,” _International journal of radiation oncology, biology, physics_ , vol. 89, pp. 191–198, 2014.
  * [46] B. Schultze, P. Karbasi, C. Sarosiek, G. Coutrakon, C. E. Ordoñez, N. T. Karonis, K. L. Duffin, V. A. Bashkirov, R. P. Johnson, K. E. Schubert, and R. W. Schulte, “Particle-tracking proton computed tomography—data acquisition, preprocessing, and preconditioning,” _IEEE Access_ , vol. 9, pp. 25 946–25 958, 2021.
  * [47] F. P. Preparata and M. I. Shamos, _Computational Geometry: An Introduction_. Springer Science & Business Media, 2012.
  * [48] M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars, _Computational Geometry: Algorithms and Applications_ , 3rd ed. Springer-Verlag, TELOS, Santa Clara, CA, USA, 2008.
  * [49] D. A. Low, P. J. Parikh, W. Lu, J. F. Dempsey, S. H. Wahab, J. P. Hubenschmidt, M. M. Nystrom, M. Handoko, and J. D. Bradley, “Novel breathing motion model for radiotherapy,” _International Journal of Radiation Oncology, Biology, Physics_ , vol. 63, pp. 921–929, 2005.
  * [50] D. A. Low, B. M. White, P. P. Lee, D. H. Thomas, S. Gaudio, S. S. Jani, X. Wu, and J. M. Lamb, “A novel CT acquisition and analysis technique for breathing motion modeling,” _Physics in Medicine and Biology_ , vol. 58, pp. L31–L36, 2013.
  * [51] M. Guo, G. Chee, D. O’Connell, S. Dhou, J. Fu, K. Singhrao, D. Ionascu, D. Ruan, P. Lee, D. A. Low, J. Zhao, and J. H. Lewis, “Reconstruction of a high-quality volumetric image and a respiratory motion model from patient CBCT projections,” _Medical Physics_ , vol. 46, pp. 3627–3639, 2019\.
  * [52] G. Chee, D. O’Connell, Y. M. Yang, K. Singhrao, D. A. Low, and J. H. Lewis, “McSART: an iterative model-based, motion-compensated SART algorithm for CBCT reconstruction,” _Physics in Medicine and Biology_ , vol. 64, p. 095013, 2019.
  * [53] J. Eom, X. G. Xu, S. De, and C. Shi, “Predictive modeling of lung motion over the entire respiratory cycle using measured pressure-volume data, 4DCT images, and finite-element analysis,” _Medical Physics_ , vol. 37, pp. 4389–4400, 2010.
  * [54] J. Ehrhardt, R. Werner, A. Schmidt-Richberg, and H. Handels, “Statistical modeling of 4D respiratory lung motion using diffeomorphic image registration,” _IEEE Transactions on Medical Imaging_ , vol. 30, pp. 251–265, 2011.
  * [55] Y. Zhang, J. Yang, L. Zhang, L. E. Court, P. A. Balter, and L. Dong, “Modeling respiratory motion for reducing motion artifacts in 4D CT images,” _Medical Physics_ , vol. 40, p. 041716 (13pp), 2013.
  * [56] M. Buhmann, _Radial Basis Functions: Theory and Implementations_. Cambridge University Press, 2003.
  * [57] R. M. Lewitt, “Alternatives to voxels for image representation in iterative reconstruction algorithms,” _Physics in Medicine and Biology_ , vol. 37, pp. 705–716, 1992.
  * [58] ——, “Multidimensional digital image representations using generalized Kaiser–Bessel window functions,” _Journal of the Optical Society of America A_ , vol. 7, pp. 1834–1846, 1990.
  * [59] G. T. Herman, “Basis Functions in Image Reconstruction From Projections: A Tutorial Introduction,” _Sensing and Imaging_ , vol. 16, pp. 1–21, 2015\.
  * [60] S. Matej and R. Lewitt, “Efficient 3D grids for image reconstruction using spherically symmetric volume elements,” _IEEE Transactions on Nuclear Science_ , vol. 42, pp. 1361–1370, 1995.
  * [61] ——, “Practical considerations for 3-D image reconstruction using spherically symmetric volume elements,” _IEEE Transactions on Medical Imaging_ , vol. 15, pp. 68–78, 1996.
  * [62] H. Heaton, C. Ordonez, Y. Censor, and R. Schulte, “Implementation of blob basis functions in proton CT reconstruction,” _Poster presented at the Joint Mathematics Meetings (JMM) of the MAA and the AMS_ , Meeting #1116, January 6-9, 2016, Seattle, WA, USA. http://math.haifa.ac.il/yair/webpage-after-200915/JMM-Blob-Poster-2016.pdf.
  * [63] Y. Benkarroum, G. T. Herman, and S. W. Rowland, “Blob parameter selection for image representation,” _Journal of the Optical Society of America A_ , vol. 32, pp. 1898–1915, 2015.
  * [64] S. Penfold, R. Zalas, M. Casiraghi, M. Brooke, Y. Censor, and R. Schulte, “Sparsity constrained split feasibility for dose-volume constraints in inverse planning of intensity-modulated photon or proton therapy,” _Physics in Medicine and Biology_ , vol. 62, pp. 3599–3618, 2017.
  * [65] S. Penfold, M. Casiraghi, T. Dou, R. Schulte, and Y. Censor, “A feasibility-seeking algorithm applied to planning of intensity modulated proton therapy: A proof of principle study,” _Medical Physics_ , vol. 42, p. 3338, 2015.
  * [66] R. Hesse, D. R. Luke, and P. Neumann, “Alternating projections and Douglas-Rachford for sparse affine feasibility,” _IEEE Transactions on Signal Processing_ , vol. 62, pp. 4868–4881, 2014.
  * [67] Y. Censor and T. Elfving, “A multiprojection algorithm using Bregman projections in a product space,” _Numerical Algorithms_ , vol. 8, pp. 221–239, 1994.
  * [68] Y. Censor, T. Elfving, N. Kopf, and T. Bortfeld, “The multiple-sets split feasibility problem and its applications for inverse problems,” _Inverse Problems_ , vol. 21, pp. 2071–2084, 2005.
  * [69] Y. Censor, A. Gibali, and S. Reich, “Algorithms for the split variational inequality problem,” _Numerical Algorithms_ , vol. 59, pp. 301–323, 2012\.
  * [70] M. Brooke, Y. Censor, and A. Gibali, “Dynamic string-averaging CQ-methods for the split feasibility problem with percentage violation constraints arising in radiation therapy treatment planning,” _International Transactions in Operational Research_ , 2020, https://onlinelibrary.wiley.com/doi/full/10.1111/itor.12929.
  * [71] C. Byrne, “Iterative oblique projection onto convex sets and the split feasibility problem,” _Inverse Problems_ , vol. 18, pp. 441–453, 2002.
  * [72] Y. Censor, “An automatic relaxation method for solving interval linear inequalities,” _Journal of Mathematical Analysis and Applications_ , vol. 106, pp. 19–25, 1985.
  * [73] K. Maass, M. Kim, and A. Aravkin, “A nonconvex optimization approach to IMRT planning with dose-volume constraint,” 2020, preprint. https://arxiv.org/pdf/1907.10712.pdf.
  * [74] V. A. Bashkirov, R. P. Johnson, H. F.-W. Sadrozinski, and R. W. Schulte, “Development of proton computed tomography detectors for applications in hadron therapy,” _Nuclear Instruments & Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_, vol. 809, pp. 120–129, 2016.
  * [75] M. Witt, B. Schultze, R. Schulte, K. Schubert, and E. Gomez, “A proton simulator for testing implementations of proton CT reconstruction algorithms on GPGPU clusters,” in _Proceedings of the IEEE Nuclear Science Symposium & Medical Imaging Conference (NSS/MIC) 2012, Anaheim, CA, USA_, 2012, pp. 4329–4334.
  * [76] B. Schultze, P. Karbasi, V. Giacometti, T. Plautz, K. E. Schubert, and R. W. Schulte, “Reconstructing highly accurate relative stopping powers in proton computed tomography,” in _Proceedings of the IEEE Nuclear Science Symposium & Medical Imaging Conference (NSS/MIC), San Diego, CA, USA, 2015_, 2015, pp. 1–3, doi:10.1109/NSSMIC.2015.7582218.
  * [77] R. P. Johnson, V. A. Bashkirov, G. Coutrakon, V. Giacometti, P. Karbasi, N. T. Karonis, C. Ordoñez, M. Pankuch, H. F.-W. Sadrozinski, K. E. Schubert, and R. W. Schulte, “Results from a prototype proton-CT head scanner,” _Physics Procedia_ , vol. 90, pp. 209–214, 2017.
  * [78] P. Karbasi, R. Cai, B. Schultze, H. Nguyen, J. Reed, P. Hall, V. Giacometti, V. Bashkirov, R. Johnson, N. Karonis, J. Olafsen, C. Ordonez, K. E. Schubert, and R. W. Schulte, “A highly accelerated parallel multi-GPU based reconstruction algorithm for generating accurate relative stopping powers,” in _2017 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), Atlanta, GA, USA_ , 2017, pp. 1–4, doi:10.1109/NSSMIC.2017.8532871.
  * [79] R. P. Johnson, V. A. Bashkirov, L. DeWitt, V. Giacometti, R. F. Hurley, P. Piersimoni, T. E. Plautz, H. F.-W. Sadrozinski, K. E. Schubert, R. W. Schulte, B. E. Schultze, and A. Zatserklyaniy, “A fast experimental scanner for proton CT: Technical performance and first experience with phantom scans,” _IEEE Transactions on Nuclear Science_ , vol. 63, pp. 52–60, 2016\.
  * [80] P. Karbasi, “Robust and efficient methods for proton computed tomography,” Ph.D. dissertation, Electrical and Computer Engineering, Baylor University, Baylor, TX, USA, 2018.
  * [81] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand, F. Behner, L. Bellagamba, J. Boudreau, L. Broglia, A. Bruneng, H. Burkhardt, S. Chauvie, J. Chuma, R. Chytracek, and G. Coopermana, “Geant4 - a simulation toolkit,” _Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 506, pp. 250–303, https://doi.org/10.1016/S0168-9002(03)01368-8.
  * [82] V. Giacometti, “Modelling and Improvement of Proton Computed Tomography,” Ph.D. dissertation, School of Physics, University of Wollongong, Australia, 2016, https://ro.uow.edu.au/theses/4903.
  * [83] J. Perl, J. Shin, J. Schümann, B. Faddegon, and H. Paganetti, “TOPAS: An innovative proton Monte Carlo platform for research and clinical applications,” _Medical Physics_ , vol. 39, pp. 6818–6837, 2012.
  * [84] P. Piersimoni, J. Ramos-Méndez, T. Geoghegan, V. A. Bashkirov, R. W. Schulte, and B. A. Faddegon, “The effect of beam purity and scanner complexity on proton CT accuracy,” _Medical Physics_ , vol. 44, pp. 284–298, 2017.
  * [85] P. Piersimoni, B. A. Faddegon, J. Ramos-Méndez, R. W. Schulte, L. Volz, and J. Seco, “Helium CT: Monte Carlo simulation results for an ideal source and detector with comparison to proton CT,” _Medical Physics_ , vol. 45, pp. 3264–3274, 2018.
  * [86] V. Giacometti, S. Guatelli, M. Bazalova-Carter, A. B. Rosenfeld, and R. W. Schulte, “Development of a high resolution voxelised head phantom for medical physics applications,” _Physica Medica_ , vol. 33, pp. 182–188, 2017\.
  * [87] S. Guatelli and S. Incerti, “Monte Carlo simulations for medical physics: From fundamental physics to cancer treatment,” _Physica Medica_ , vol. 33, pp. 179–181, 2017.
  * [88] V. Giacometti, V. A. Bashkirov, P. Piersimoni, S. Guatelli, T. E. Plautz, H. F.-W. Sadrozinski, R. P. Johnson, A. Zatserklyaniy, T. Tessonnier, K. Parodi, A. B. Rosenfeld, and R. W. Schulte, “Software platform for simulation of a prototype proton CT scanner,” _Medical Physics_ , vol. 44, pp. 1002–1016, 2017.
  * [89] R. Schulte, K. Schubert, and Y. Censor, “Computational tools for real-time adaptive particle therapy,” 2013, this is a research grant from the United States-Israel Binational Science Foundation (BSF), Grant Number 2013003, 10/2014-9/2019.
  * [90] R. Schulte and Y. Censor, “Intensity modulation and proton computed tomography in proton radiation therapy planning,” 2009, this is a research grant from the United States-Israel Binational Science Foundation (BSF), Grant Number 2009012, 10/2010-9/2013.

|  Yair Censor received his B.Sc., M.Sc. and D.Sc. degrees in mathematics from
the Technion-Israel Institute of Technology in Haifa in 1967, 1969 and 1975,
respectively. From 1977 to 1979, he was Research Assistant Professor of
Computer Science at the State University of New York (SUNY) at Buffalo, NY,
USA, and Research Associate with the Medical Image Processing Group (MIPG)
there. In 1979, Dr. Censor joined the Department of Mathematics at the
University of Haifa in Israel, where he has been Full Professor since 1989 and
became Professor Emeritus in 2012. He was a visiting professor at the
University of Pennsylvania in Philadelphia; in Linköing University in Sweden;
in the Instituto de Matematica Pura e Aplicada in Rio de Janeiro, Brazil; in
the Beijing International Center of Mathematics at the Peking University in
Beijing, China; and in the Graduate Center of the City University of New York.
Dr. Censor works in computational mathematics, where his interests include
optimization theory, inverse problems, optimization techniques in image
reconstruction from projections, and in intensity-modulated radiation therapy
(IMRT). Dr. Censor has published over 160 research articles in refereed
scientific journals, conference proceedings and as book chapters. He co-
authored with S.A. Zenios the book: Parallel Optimization: Theory, Algorithms,
and Applications, Oxford University Press, New York, 1997. For this book he
received, together with Professor Zenios, the 1999 ICS (INFORMS Computing
Society) Prize for research excellence in the interface between operations
research and computer science. Dr. Censor served as Associate Editor for the
journal IEEE Transactions on Medical Imaging from 1992 to 2004.
---|---
|  Keith E. Schubert received his B.S. degree in general engineering from the
University of Redlands in 1991, his M.S. in electrical engineering from UCLA
in 1992. After working for a few years, he earned his Ph.D. in electrical and
computer engineering from UCSB in 2003 for work on computational methods of
uncertain systems and robustness in estimation and control. From 1992-1996 he
worked for Northrop-Grumman as a design engineer, and since 1996 he has worked
as a consulting engineer for various companies. In 2000, he took a visiting
faculty position at the University of Redlands in mathematics and computer
science for two years, then was hired as an assistant professor in 2002 at CSU
San Bernardino, where he became an associate professor in 2006, then a full
professor in 2010. At CSUSB he was the lead on starting their first
engineering program and was part of the team that started the Bioinformatics
program. In 2013, he joined Baylor University in Waco, Tx, where he is
currently a full professor in electrical and computer engineering. His
research is in devices and computational methods for robust imaging
particularly for biomedical applications. Dr. Schubert has several patents and
served as Principal Investigator on the NIH R01 grant to develop proton CT for
biomedical applications from 2011-2016 and as part of the research team for a
planning grant (P20) for the development of ion therapy research in the United
States from 2015-2017. In 2013, Dr. Schubert was one of the team leaders of a
National Geographic research expedition to Cueva de Villa Luz in Mexico, where
he studied patterned biological grows in the sulfuric acid caves, which was
featured in the July 2014 cover article of National Geographic Magazine.
---|---
|  Reinhard W. Schulte (M’12) received his diploma degree (M.S. equivalent) in
physics (Diplom Physiker, Dipl. Phys.) from the University of Dortmund,
Germany, in 1978, and his doctorate in medicine (Doctor of Medicine, Dr. med.)
from the University of Cologne, Germany, in 1986. He completed his residency
in radiology and radiation oncology at Hannover Medical School in 1989. He was
a postdoctoral fellow at the first hospital-based proton therapy center at
Loma Linda University Medical Center (LLUMC) in California from 1990-1994. In
1994 he accepted a clinical research position at LLUMC. In 1996, he became an
Assistant Clinical Professor in the School of Medicine and the Department of
Radiation Medicine at LLUMC, In 2007, he was promoted to Associate Professor,
and since 2013, he has been a Full Professor at Loma Linda University, now in
the Department of Basic Sciences, Division of Biomedical Engineering Sciences.
He has been successful in developing novel instrumentation and algorithmic
approaches for proton therapy. His current research interests include the
clinical development of proton CT and nanodosimetry to improve proton and ion
therapy. He is an Associate Editor of the journal Medical Physics and holds
several patents. Prof. Schulte was awarded a habilitation fellowship by the
German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) from
1989-1991. He is a member of the American Society of Radiation Oncology
(ASTRO), the European Society of Therapeutic Radiation Oncology (ESTRO), the
American Association of Physicists in Medicine (AAPM), the German Society for
Boron Neutron Capture Therapy, and the Radiation Research Society. He was
awarded the ESTRO Calergo Award in 1990. Prof. Schulte has served as Principal
Investigator on several NIH grants, including an R01 grant to develop proton
CT for biomedical applications from 2011-2016 and a planning grant (P20) for
the development of ion therapy research in the United States from 2015-2017.
---|---