# The Colibri Telescope Array for KBO Detection through Serendipitous Stellar
Occultations: a Technical Description

Michael J. Mazur 1,∗, Stanimir Metchev 1,2 Rachel A. Brown 1,2, Ridhee Gupta
3, Richard Bloch 4, Tristan Mills 1, Emily Pass 5

###### Abstract

## 1

We present the technical design, construction and testing of the Colibri
telescope array at Elginfield Observatory near London, Ontario, Canada. Three
50-cm telescopes are arranged in a triangular array and are separated by
110-160 metres. During operation, they will monitor field stars at the
intersections of the ecliptic and galactic plane for serendipitous stellar
occultations (SSOs) by trans-Neptunian objects (TNOs). At a frame rate of 40
frames per second (fps), Fresnel diffraction in the occultation light curve
can be resolved and, with coincident detections, be used to estimate basic
properties of the occulting object. Using off-the-shelf components, the
Colibri system streams imagery to disk at a rate of 1.5 GB/s for next-day
processing by a custom occultation detection pipeline.

The imaging system has been tested and is found to perform well, given the
moderate site conditions. Limiting magnitudes at 40 fps are found to be about
12.1 (temporal SNR=5, visible light Gaia $G$ band) with time-series standard
deviations ranging from about 0.035 mag to >0.2 mag. SNR is observed to
decrease linearly with magnitude for stars fainter than about $G=9.5$ mag.
Brighter than this limit, SNR is constant, suggesting that atmospheric
scintillation is the dominant noise source. Astrometric solutions show errors
typically less than $\pm$0.3 pixels (0.8 arc seconds) without a need for high-
order corrections.

##  2 Keywords:

robotic telescopes, Kuiper Belt Objects, occultation, photometry, Fresnel
diffraction

##  3 Introduction

###  3.1 The Undiscovered Population of Kilometre-Sized Outer Solar System
Objects

Beyond Neptune orbits a population of up to $10^{11}$ objects with sizes of a
kilometre or larger (Roques and Moncuquet, 2000), known as trans-Neptunian
objects (TNOs). It can be divided into three main subgroups: Kuiper Belt
objects (KBOs), scattered disk objects, and the detached objects such as the
Sednoids and the distant Oort Cloud objects (Delsanti and Jewitt, 2006).
Extending from about 30 AU to about 50 AU from the sun, the Kuiper Belt
contains objects with predominantly low inclinations that are dynamically
stable compared to the broader group of TNOs. The KBOs can be subdivided into
two groups known as the dynamically cold and hot populations. The cold
population typically have low inclination (<5°) orbits and are believed to
represent a dynamically primordial population (Levison and Stern, 2001). The
cold population is also expected to have the highest sky surface density,
since it most closely follows the solar system ecliptic plane.

Compared to asteroids, for which the size-frequency distribution is reasonably
well understood, less is known about this population of distant solar system
objects. Because of their distance from Earth, smaller KBOs can be difficult
or impossible to observe directly. Of the more than one thousand KBOs
discovered by direct imaging, only a small fraction (3%) are less than 25 km
in size and none are less than 7 km. As a result, the size-frequency
distribution of KBOs and, more generally TNOs, is poorly defined for smaller
objects.

###  3.2 Indirect Detection through Serendipitous Stellar Occultations

Using stellar occultations as a method for the indirect detection of
‘invisible’ bodies in the Solar System was suggested by Bailey (1976). Using a
1 m telescope with a vidicon imager and a visual ($V$-band: $\approx$505–595
nm) limiting magnitude of 16, Bailey suggested that a 1000-star field would
have an occultation event every 11 hours. His analysis, however, is optimistic
as it does not consider the effects of line-of-sight Fresnel diffraction nor
fully explore the temporal limit imposed by a 10 Hz sample rate.

A theoretical description of Fresnel diffraction during stellar occultations
by small bodies was discussed by Roques et al. (1987). Fresnel diffraction is
used to describe near-field diffraction and, given the correct geometries and
source/object sizes, can be used to describe the observed light curves from
KBO occultation events. Roques et al. (1987) examined both the theory and
implementation of diffraction models and forms the basis of subsequent
observational studies in the field. Roques and Moncuquet (2000) further
explored the possibility of detecting small bodies in the Solar System. Their
work focused on sub-kilometre KBOs and predicted the possibility of a few to
several tens of occultation events per night with 2-m to 8-m class telescopes
at a visual limiting magnitude of 10 mag $\leq V\leq$ 12 mag. With today’s
instrumentation, however, similar limiting magnitudes are achievable with
smaller telescopes operating at higher frame rates.

In order to sufficiently resolve a KBO occultation light curve, we need to
image at relatively high frame rates. For a 1 km KBO at 40 AU observed at
solar opposition, its speed relative to the observer is about 25 km/s. If this
object occults a star with a 1 km projected diameter, the occultation will
have a duration $\Delta T\approx 80$ ms if we consider the event purely
geometrically. However, the diffractive broadening of the geometric shadow
prolongs the event by a factor of 2–3. The effect is more pronounced for
smaller stellar disks and better photometric precision (Roques et al., 1987;
Roques and Moncuquet, 2000).

Taking into account relative velocities, Fresnel geometry, and telescope
sensitivity, Bickerton et al. (2009) find that the optimal sampling rate for
detecting serendipitous stellar occultations by KBOs at visible wavelengths is
40 fps, with observations toward solar opposition. These are the parameters
that we use to set the technical requirements for the Colibri telescope array.
These parameters are preliminary. Until reliable KBO occultation detections
become routine, the optimal trade-off between telescope size, sampling rate,
observing wavelength, and observing geometry has yet to be empirically
validated.

Indeed, the Bickerton et al. (2009) study details a much longer set of
assumptions to project KBO occultation detection rates. We defer a discussion
of the expected event rate specific to the Colibri experiment to a future
publication (Metchev et al. 2022, in preparation). In the meantime, we note
that some recent and on-going experiments (Section 3.3) validate our choice of
telescope system and operating mode.

###  3.3 Other Previous or Planned Experiments

Several prior experiments have reported serendipitous stellar occultations by
KBOs. These have often been based on sub-optimal data sets, in some cases
acquired for a different science goal. Chang et al. (2006) analyse 90 hours of
archival x-ray monitoring observations of a single source, Cygnus X-1, with
the Rossi X-ray Timing Explorer (RXTE) satellite at 500 fps. The 58 candidate
occultations reported in their data are by far the largest in a single data
set. However, subsequent re-analyses of the data and their statistical
interpretation by Jones et al. (2006) and Bickerton et al. (2008) have put
these detections into doubt. Schlichting et al. (2009, 2012) report two
different candidate occultations from visible-light guiding operations at 40
fps from over 20 years of observatinos with the Hubble Space Telescope. These
events do bear the hallmarks of the expected Fresnel diffraction pattern of
stellar occultations by kilometre-sized KBOs (e.g., Figure 1 of Schlichting et
al., 2009). Detections of similar events with other facilities would confirm
them as representative of this phenomenon.

Early observations designed specifically for the detection of serendipitous
stellar occultations have also yielded some candidate detections and mixed
results. Roques et al. (2006) report three candidate events in 11 h of dual-
band visible wavelength monitoring of two stars at 45 fps with the 4.2 m
William Herschel Telescope. Bickerton et al. (2008) discuss that while the
rate of events in this study is in line with expectations, their statistical
significance is low. More recently, Arimatsu et al. (2019) report a single
candidate event using a pair of 28 cm amateur optical telescopes, in a 60 h
observation at 15.4 fps in the course of the Organized Autotelescopes for
Serendipitous Event Survey (OASES; Arimatsu et al., 2017). The sub-optimal
(<40 fps) cadence of the observations and the low (<10) SNR of the four
individual measurements that constitute the candidate event leave enough room
for it to be a false positive. Nonetheless, the OASES setup and its use of
commercial hardware and a rapid-imaging complementary
metal–oxide–semiconductor (CMOS) camera are promising for designing large-
scale serendipitous stellar occultation surveys.

Most significantly, the Taiwanese-American Occultation Survey (TAOS; Lehner et
al., 2009) was specifically designed to identify $\sim$1 km-diameter objects
beyond the orbit of Neptune, and to measure the size distribution of KBOs with
diameters between 0.5–30 km. Seven years of visible-light monitoring with
initially three, and then four 50 cm telescopes with TAOS yielded no
occultation detections (Zhang et al., 2013). This was attributed to a lower-
than-expected event rate, and also the relatively slow (5 fps) sampling of the
cameras. A follow-up experiment, TAOS II (Lehner et al., 2012), is in the
final stages of development, and will use three 1.3 m telescopes imaging at 20
fps. Much like OASES, and the herein described Colibri Telescope Array, TAOS
II will use a CMOS-type visible-light camera.

##  4 Colibri Hardware

Since stellar occultations by KBOs tend to be short-lived, they are difficult
to observe. Any experiment to identify these transient events needs to not
only rapidly (at 40 fps) image background stars but should also provide a
mechanism for confirmation. The three telescopes of the Colibri array are set
up in a triangular pattern with spacings of between 110 m and 160 m from each
other (Figure 1). This arrangement allows us to rule out the possibility of
atmospheric twinkling for coincident events but is not sufficient for
size/distance determination by differential transit timing. A full description
of the technical specifications of the Colibri array is given in Table 1.

Figure 1: A satellite image showing the location of the three Colibri telescopes (stars) and the control center (CC) located in the basement of a house on the site. Maps Data: Google. Imagery ©2022 First Base Solutions, Maxar Technologies. Table 1: Description of the three-telescope Colibri array.  | Specification
---|---
Location 1 (Lat., Long.) | 43.193365°N, 81.316090°W
Location 2 (Lat., Long.) | 43.192415°N, 81.316467°W
Location 3 (Lat., Long.) | 43.193108°N, 81.318033°W
Telescope baselines | 110–160 m
Telescopes | Hercules 50 cm f/3 w/ Wynne Correctors
Cameras | FLI Kepler KL4040
Sensors | GSense 4040
Pixels | 4096x4096, 9$\mu$m pitch
Field-of-View | 1.43°x1.43°
Pixel Scale ($2\times 2$ binning) | 2.52′′ per binned pixel
Digitization | low/high gain @ 12-bit
Readout Noise | 3.7 $\mathrm{e^{-}}$
Dark Current @ -20°C | <0.15 $\mathrm{e^{-}}$/pixel/s
Typical Gain | low= 19 $\mathrm{e^{-}}$/ADU, high= 0.8 $\mathrm{e^{-}}$ /ADU
Quantum Efficiency | >70$\%$ (460 nm – 680 nm)
Framerate @ $2\times 2$ binning | >40 fps
Mag. Limit (SNR=3) in 25 ms | $G=12.5$
Mag. Limit (SNR=5) in 25 ms | $G=12.1$

###  4.1 Telescopes and Domes

The light-weight carbon fiber optical tube assemblies support 50 cm f/3
cellular mirrors cast from Schott Boro 33 glass. The telescope manufacturer is
Hercules Telescopes, which both built the OTAs and cast the mirrors. The
mirrors were then ground and figured by an external contractor. ASA Wynne
correctors correct field aberrations for imaging at prime focus. FLI Kepler
KL4040 cameras are attached to a focus assembly that is controlled by a
Seletek Platypus controller. The controller is connected to our local VLAN
where it is accessible to the control computers. Absolute timing is achieved
by a Garmin GPS attached to each camera.

The telescopes are attached to Astro-Physics AP1600 GTO equatorial mounts that
sit on custom-built steel piers (Figure 2). The steel piers are, in turn,
attached to concrete piers extending to about 6 metres below the dome floor.
Protection from the elements is provided by ExploraDome EDII 12′ plastic domes
that have a number of modifications to increase their reliability and allow
them to be operated robotically, under the icy and windy conditions common at
the site. The original dome control electronics proved to be unreliable and so
were replaced with MaxDome II controllers. The dome controllers are connected
to USB-over-IP appliances for dome control over much greater distances than
would otherwise be possible.

Figure 2: The Hercules 50-cm telescope is mounted on an AP1600 GTO mount
inside of a 12′ automated dome.

Environmental conditions are monitored on-site by a number of different
systems. A central Davis Pro weather station provides an overview of the
weather while each dome will be equipped with its own networked environmental
station that measures inside/outside temperatures and humidities as well as
ambient light levels. Additionally an Arduino-based IR cloud monitor and a
dedicated camera continuously aimed at the star Polaris give feedback on
seeing conditions. Python scripts talk to each of the devices, parsing and
collating the data into a format that can be read by the telescope control
software.

###  4.2 Imaging Cameras

Traditional charge-coupled device (CCD) cameras – available in large formats
at reasonable costs – are still a common option for many astronomical imaging
applications. Although most off-the-shelf solutions use high-speed USB
interfaces, their frame rates are typically no faster than a few frames per
second. As this is still too low for resolving most KBO occultations,
commercial CCD cameras do not meet our requirements.

There are other cameras, using electron-multiplying CCDs (EMCCDs), that often
have higher frame rates and better noise characteristics than regular CCD
cameras. By using an electron-multiplying stage prior to their output
amplifier, EMCCD cameras can achieve gains of more than 1000 while maintaining
a read noise on the order of a few electrons. Their cost, however, tends to be
high and their sensor sizes small. A high-sensitivity EMCCD camera borrowed
from the Western Meteor Physics Group was tested, but its high cost, moderate
frame rate (17 fps), and small chip dimensions (13 mm x 13 mm) effectively
ruled it out as a viable option. The testing did, however, allow the Colibri
processing pipeline to be validated on real data and reinforce the need for
higher frame rates and a large field of view (Pass et al., 2018).

Over the past few years, there has been a shift within the imaging community
from CCD sensors to CMOS sensors. While CMOS-based cameras have been found in
cameras and mobile phones for years, scientific CMOS (sCMOS) cameras are
relatively new to the astronomical market. Their architecture differs from
CCDs as, instead of an output register through which pixel data passes
serially, they have an output amplifier for every column so that data from
individual pixels can be read out in parallel. This has the advantage of
significant speed increases over CCDs, but typically comes with the tradeoff
of having a rolling shutter instead of a global shutter. Although global
shuttering has the desired characteristic of instantaneous capture of an
entire frame, a rolling shutter is only a liability when objects are moving
across the field of view. Compensations for timing offsets between rows can
easily be made during processing if deemed necessary. With a 25 ms exposure
time and 2048 (2x2 binning) rows being read out, the row-to-row timing offset
is 1.2 $\mu$s. Similar to the accuracy of our time stamp, it is not
significant as long as the pointing differences between the three telescopes
is small. Considering all of this, sCMOS-based cameras were the only
technology that could meet the speed, field-of-view, and cost requirements of
the Colibri project.

Although sCMOS-based cameras were the preferred choice, no suitable cameras
were commercially available at the time of project conception. The cameras
that were on the market were either equipped with small chips or delivered
frame rates much lower than the required 40 fps. At the time, FLI was
developing their Kepler series cameras which could be equipped with a large-
format sCMOS sensor and, most importantly, would be able to image at frame
rates >30 fps over a fiber connection. Because of this, FLI was chosen to
supply three KL4040 cameras for the project.

The cameras have a number of possible modes of operation. They can be run in
image capture mode where data is buffered and written to disk in several
different formats, or they can stream straight from camera to disk. When
running in the first mode, the cameras can be connected to the PC via either a
USB3.0 or a fiber connection. When running over USB, speed is limited and we
find a maximum frame rate of about 15 fps when 2x2 binned. A fiber connection
allows for much higher frame rates (40+ fps), however the raw-to-fits format
conversion when writing to disk slows things considerably after the camera’s
buffer has been filled. To allow for the required continuous high-speed
imaging, data from the cameras is written directly to a RAID from a raw data
stream with no format conversion. The raw format of the stream is one where
each image container holds a header and both the low- and high-gain exposure
data. The header takes up the first 246 bytes and contains necessary
information about the exposure such as integration time and GPS time. After
the header, the 12-bit high- and low-gain image data is interleaved by rows,
giving a resultant raw image size of 12.6 MB. The two images can then be split
and combined to create a single high dynamic range (HDR) 16-bit image. As the
creation of an HDR image adds significant processing time, we currently
extract only the high-gain image for processing. Because the use of the high-
gain image gives roughly 2 magnitudes greater depth compared to the low-gain
image without saturating stars in our target fields, trading off dynamic range
and a slight increase in noise for processing speed has only a minor effect on
our ability to detect SSOs with our pipeline.

###  4.3 Network and Control Computers

Two fiber networks connect the imaging equipment in the domes with control
computers in a house located 270 m (along the fiber) from the furthest dome.
The first network uses 1000BASE-SX transceivers over single-mode fiber and is
used for communications with the telescope, focuser, guide and security
cameras, and dome controller. The second network is used exclusively by the
imaging camera. To keep up with the camera’s high frame rate, variable bit-
rate 40GBASE-CSR QSFP+ transceivers run at speeds of 8 Gb per second on each
of four channels over 8 strands of OM3 multi-mode fiber. Speed tests have
shown that frame rates up to 50 fps with 2x2 binned images are achievable, but
are at the upper limit of our four-disk RAIDs.

Each dome has its own dedicated control computer to manage equipment
operations, scheduling, and imaging tasks. The PCs are running Windows 10 with
observatory control managed by ACP software. ACP allows us to connect to any
device with ASCOM drivers which makes integration of off-the-shelf equipment
relatively straight-forward. Technically, the PCs are mid-range machines with
3.2 GHz i7-8700 processors and 32 GB of RAM. To allow for image streaming to
disk at rates of up to 4 Gb per second, each computer has a striped RAID of 4
10 TB Seagate Exos hard disks. This gives a peak performance of about 5 Gb per
second with a total storage capacity of 36 TB. At 40 fps this gives enough
space for about 20 hours – roughly 2 to 3 nights – of data.

##  5 Colibri Data Handling

###  5.1 Data Acquisition

Since the streaming mode of camera operation is not available through the FLI
ASCOM driver, we have written our own camera control code in C#. This
lightweight application runs from the command line and gives the user access
to the main camera settings while running in streaming mode. Command line
usability also means that the program is accessible from within ACP’s
scripting environment.

The lack of ASCOM camera support for our operations also means that we need an
alternative to ACP’s scheduler. We accomplished this by writing JavaScript
code to handle rudimentary scheduling that could be run from ACP. After
initiation, the script checks for the presence of weather data. If present,
the script continues. Otherwise, the user is given the opportunity to either
quit or continue without the safety of the weather station data. The next step
creates the data directories for that night’s data and collects a set of 50
bias images. The script then calculates sunset and sunrise times, as well as
the position of the moon relative to a pre-determined list of stellar fields
of observation. Details on the assembly of this pre-determined set of stellar
fields are deferred to a future publication (Metchev et al. 2022, in
preparation).

We iterate through each field in the list and estimate the number of visible
stars according to airmass in 6 minute steps from sunset to sunrise. The
fields are then filtered based on elevation (must be >10°) and moon angle
(must be >15°), with ones that do not satisfy these criteria removed from the
list of fields for the night. The remaining fields are then ranked, according
to the estimated number of stars visible, for each 6 minute window. Once this
has been done at the beginning of a night, a sunset check is made. If the sun
is less than 12° below the horizon, the system pauses until the sun has passed
this limit.

Once it has been determined that it is dark enough to open, a weather-check is
made and, if it is safe, the dome is opened, homed, and synced to the
telescope. After another set of biases are collected, the telescope then slews
to the first field and begins capturing data in subsets of 2400 images. Once
another field surpasses the current field, the telescope is slewed to the new
location and imaging begins once again. This is repeated until either sunrise
or a weather alert shuts the system down.

###  5.2 Preliminary Data Processing

After being written to disk, the images can be viewed individually with either
FLI’s Pilot software or with custom python scripts for viewing and conversion
to FITS format. Batch conversion to FITS format averages about 20 ms per
image. This can be done after a night of observations, but effectively takes
nearly as long as the observations themselves. It is therefore an important
limitation as we aim to process data during daylight hours. Our data storage
allows up to a few consecutive clear nights of observations, which is
generally sufficient given the $\sim$30% fraction of clear nights at
Elginfield Observatory.

As with CCD imagers, sCMOS cameras should be corrected for bias, flat-field,
and dark current prior to quantitative/absolute work being done. Although the
noise specifications are quite low, the sensors in our cameras have a
considerable amount of bias structure that is especially visible during rapid
imaging. As such, it must be removed by subtraction of a median combined set
of bias frames.

With the short (25 ms) exposures that are typical for Colibri operation, the
low dark current (0.15 $\mathrm{e^{-}/pixel/s}$) of the sensor means that we
can ignore separate dark current removal. Even then, because of the way that
we collect the biases (25 ms, shutter closed), dark current is effectively
removed during bias subtraction. Collection of a traditional bias with a 0 ms
exposure is not possible currently due to limitations with the camera
hardware.

The fast, wide-field nature of our optical system lends itself to vignetting
on the image plane. This is visible as darkening outwards from the centre of
the field and is easily corrected with a proper flat-field image. For stellar
extraction and the relative photometry that we perform in our workflow,
however, flat-fielding is not strictly necessary and is not applied during the
first stage of processing due to time constraints.

###  5.3 Occultation Detection Pipeline

Although the detection pipeline (Figure 3) follows the general description
given by Pass et al. (2018), changes to the equipment specifications have
necessitated modifications to the pipeline.

Figure 3: Processing pipeline flowchart showing the steps in the Colibri data
reduction pipeline.

After the initial preparation of the data has been completed, stellar sources
are extracted by running the Source Extraction and Photometry (SEP) module for
Python (Bertin and Arnouts, 1996; Barbary, 2016) on a median combination of
the first nine images from each 2400 image subset collected over the night.
The output from the SEP $extract$ function gives the coordinates for each star
with a flux at least 4 $\sigma$ above the background at the start of each
subset. After this, the average stellar drift within each subset is
calculated. Every image in a subset is then interrogated using the
$sum\\_circle$ function of SEP to produce a time series of fluxes at each
drift-corrected x,y location. The light curves are then passed to the first-
stage dip detection where stars are filtered based on their SNR (must be
greater than 5) and whether they drift out of the field. We then look for
geometric dips greater than 40% of the normalized light curve. A dip much less
than 40% is considered a common event (due to noise) and would represent a
false detection. The geometric dip detections are automatically saved as
candidates. In the case of no geometric dip, the light curve is passed to the
KBO dip detection function which convolves a Ricker wavelet designed to match
the characteristic width of a KBO occultation event with each time series. The
size of the wavelet is chosen based on the expected occultation duration given
the pointing of the telescope. This means that, with expected durations of
between 160 ms and 500 ms, Ricker wavelets of between 4 and 12 frames length
are most appropriate. We then look for the most significant dimming events in
each time series. If the dips of these events are greater than 3.75$\sigma$ of
the mean of the convolved time series, the event is passed to the next stage
of processing where the events are compared to a set of pre-calculated kernels
that model diffraction patterns for a range of different physical/dynamical
properties. Finally, successful matches from this step are checked for
correlation in time with the other two telescopes of the array and saved
offsite.

##  6 System Performance

###  6.1 Optical

When pairing a fast (f/3) mirror design with an off-the-shelf (ASA Wynne)
field corrector, the final results can be difficult to predict. The
prescription for the ASA corrector is unknown, but the manufacturer claims
that it has been designed with parabolic primaries from f/3–f/5 in mind.
Extrapolating from their sales literature, we would expect spots ranging from
about 4–12 $\mu$m (400–700 nm) on the optical axis to 12+ $\mu$m in the
corners of our sensor. Although much better than an uncorrected system, our 9
$\mu$m pixels mean that we should be able to observe the difference in spot
sizes from centre to sensor corners (Figure 4). Raytracing the as-built mirror
specs with a generic Wynne corrector prescription, however, suggests that the
design has excellent optical performance over the entire field for our
observing conditions and our 2.5′′ ($2\times 2$ binned) pixels. With seeing
often more than several arcseconds at the Elginfield site, we can tolerate
slightly defocused images while still maintaining critically sampled PSFs
across our $2\times 2$ binned images.

Figure 4: The four corners of this 10 s, $2\times 2$ binned image (2.51′′ per
binned pixel) show spot sizes visually similar to the central portion of the
image. The outer fields are 30′ square while the central square is 20′ square.

Although the coma-corrected field-of-view is claimed to be 60 mm, there is
vignetting at the field edges. This is visible in flat field imagery and
amounts to a loss of 1 magnitude in the corners of the images. As a result, a
flat-field correction should be done if we want to perform absolute
photometry. Because of Colibri’s limited processing power and the fact that we
are doing same-star differential photometry, however, we do not flat-field
correct our images within our processing pipeline.

###  6.2 Mechanical

After coarse polar alignment, the mounts were polar-aligned more precisely
using the drift-alignment technique. Pointing models to 10 degrees altitude
were then created using MaxPoint software from Diffraction Limited to quantify
mount and telescope alignment. Mount alignment was then adjusted and new
pointing models created in an iterative fashion until azimuth and altitude
misalignments were below 10′′. After alignment and pointing model corrections,
we find that the pointing performance of all telescopes is typically within a
few 10s of arcseconds of the desired position. With our large fields, this is
sufficient for ensuring that all three telescopes are monitoring the same
stars.

###  6.3 Photometry

When imaging rapidly with the KL4040 camera, background values can be seen to
fluctuate by 10s of counts. To try to understand the effect that this has on
rapid-imaging performance, we image a star-rich region of the sky (R.A.:
4.75h, Dec.: 72.75°) at our target frame rate of 40 fps without a filter. The
images are corrected for bias and then stars are extracted using SExtractor
(Bertin and Arnouts, 1996). A world coordinate system solution is obtained by
processing an initial image of the chosen stellar field through the
astrometry.net algorithm (Lang et al., 2010). Stellar positions are cross-
correlated with their Gaia Early Data Release 3 (EDR3; Gaia Collaboration et
al., 2016, 2021) coordinates. Instrumental magnitudes are then compared with
Gaia $G$-band (400–860 nm) magnitudes to derive a first-order transformation
equation (Figure 5). Using a Random Sample Consensus (RANSAC) algorithm to
estimate the best fit to our data we find very similar solutions for each of
the sigmas tested. One of the main benefits of the RANSAC solution is that
misidentified stars or those with bad instrumental magnitudes will be
classified as outliers and, as a result, have no impact on the final solution.

Stellar limiting magnitude has been defined Harris (1990) as the magnitude at
which only 50% of the objects of that magnitude are detected. By injecting
stellar sources with a range of magnitudes into real imagery, Harris counts
the number of found sources with the IRAF/DAOphot software. The method we
describe is slightly different, yet is still effective for our purposes.
Instead of injecting synthetic stars into our images, we simply compare our
extracted stars to known sources from the Gaia EDR3 catalogue to see at what
limit we stop detecting objects. Figure 5 shows the instrumental magnitude of
detected objects to a limit of 3 sigma above the background plotted against
the Gaia $G$-band magnitude of the star closest to the extracted position.
Beyond about $G=12.5$ mag, correlation between detections and real stars is
lost and we set that as our practical limiting magnitude for a 25 ms exposure.

Figure 5: Gaia $G$-band vs. instrumental magnitude of 129 sources extracted
with a detection threshold $\geq$3 sigma (left panel) and 69 sources extracted
with a detection threshold of $\geq$5 sigma (right panel)

from a 25 ms exposure.

Although an SNR of 3 is considered to be the absolute limit for most systems,
the actual limit for observing occultation events is somewhat higher. For the
detection of KBOs, we set our stellar detection threshold to 4 sigma but
realistically expect to be able to detect SSOs down to about 5 sigma above the
background. We use the fitted relation in Figure 5 to calibrate our conversion
from Colibri instrumental to Gaia $G$-band magnitudes. Since Colibri’s is a
filter-less system, we do not seek to accurately transform the data to a
standard photometric system. Our conversion is likely accurate to $\pm$0.1 mag
for most stars.

Using the same 40 fps dataset as above, we also examine the relative
photometric temporal stability of the system. To do this, we look at the
frame-to-frame residual changes in measured magnitude (instrumental or
transformed) for all stars in the field. Figure 6 shows the measurements for a
randomly-selected $G=10.1$ mag star over 1000 consecutive 25 ms exposures (25
sec sequence). The standard deviation is 0.057 mag. We use a derivative
measure, the ratio of the mean flux to the r.m.s. scatter of the flux over a
60-sec interval, as the “temporal SNR” of our flux measurements. The temporal
SNR is approximately diagnostic of the minimum detectable depth for a stellar
occultation. For example, at a temporal SNR of 10 the rms scatter of the light
curve is $1/10=1\%$, so detectable occultation would be at least $\sim$30–40%
deep. Our first dip-detection pass in our occultation detection pipeline does
indeed set a 40% threshold for detection (Section 5.3).

Figure 6: A magnitude vs. time light curve for a $G=10.1$ mag star over a 25
second window in one thousand 25 ms exposures.

Figure 7 presents the temporal SNR of objects (airmass $\approx$ 1.37) with 25
ms exposures as a function of stellar $G$ magnitude. The temporal SNR follows
a linear trend for stellar magnitudes fainter than about $G=9.5$ mag. Brighter
objects, however, are best described by a horizontal line at an SNR of about
26. This is likely an indicator that atmospheric scintillation becomes the
dominant source of noise for bright objects.

When the standard deviations of the measurements for all stars against their
mean magnitudes are plotted (Figure 7), they show a monotonically increasing
standard deviation with increasing stellar magnitude. This provides a useful
estimate of the threshold for the depth of a detectable dip at a given
magnitude. The 0.035 mag error ($G\lesssim 9.5$ mag stars), for example,
corresponds to 3.5% of the flux while a 0.2 mag error ($G=11.5$ mag stars
corresponds to 17%) of the total flux. Therefore, we can conclude that most
>20%-deep occultations of brighter stars should be visible. Occultations of
fainter stars are more likely to be detected under ideal conditions or if
$\gtrsim$50% depths occur under favourable occultation geometry. As already
detailed in Pass et al. (2018), stellar brightness will be critical for the
precision to which the parameters of a stellar occultation can be determined.

Figure 7: The temporal SNR of $G\lesssim 10$ mag objects, extracted from a
field centred at an airmass of 1.37, is roughly SNR=26. At fainter magnitudes,
the temporal SNR decreases linearly to SNR = 5 at $G\approx 12.5$ mag.

###  6.4 Astrometry

Using the SkyFit2 tool (Vida et al., 2021), a plate solution is computed using
a third-order polynomial with a first-order radial term (Figure 8). When
compared to Gaia EDR3, the results show the standard deviation of the
residuals ranging from about 0.2 to 0.25 pixels (0.5′′ to 0.63′′) with no
obvious dependence on distance from the field centre. Higher order fits have
been tested, but there is no obvious benefit to using them given the low
amount of distortion (Figure 8, upper panel) in the system and the fact that
high astrometric precision is not a strict requirement for the detection of
KBOs by stellar occultation. Although this represents a solution from a single
pointing direction (altitude=45°, azimuth=0°), we don’t find a large changes
in the plate solution at different pointings.

Figure 8: Distortion (scaled 20x) as computed using a third-order polynomial
with a first-order radial term is shown in the upper panel. The lower panel
shows the astrometric residuals after correction.

##  7 Summary

The Colibri telescope array at Elginfield Observatory in Ontario, Canada, has
been built for the purpose of detecting kilometre-sized KBOs by the
serendipitous stellar occultation method. The construction used cost-
effective, off-the-shelf components to meet the design goals of the Colibri
project.

Using sCMOS cameras, Colibri has been designed to continuously monitor the
night sky for serendipitous stellar occultation events. At 40 fps, each camera
streams imagery at 4 Gb/s over a 40G fibre link to its own dedicated RAID.
This generates up 20 TB of data per camera per night to be processed the
following day by the processing pipeline.

The photometric and astrometric performance of the system has been measured.
Photometrically, the system performs well with a limiting broad-band $G$
magnitude of about 12.5 (temporal SNR=5) in a 25 ms exposure. The relationship
between instrumental and Gaia $G$ magnitudes is linear with a 1:1 slope while
light curve analysis shows stability in the measurements within the timescale
of the 1-minute data subsets. The standard deviation in the magnitude
measurements increases with increasing magnitude, ranging from about 0.035 mag
on $G\lesssim 9.5$ mag stars to $\sim$0.2 mag at the SNR=5 limit. As this
corresponds to a roughly 3.5% to 17% variation in flux, occultations with
$\sim$50% flux depths over several consecutive exposures should be detectable
at even the faintest magnitudes. Astrometrically, the optics show some
distortion and better than $\pm$0.25 pixel ($\pm$0.6′′) errors. As a result, a
third-order polynomial plate solution is sufficient for this project.

The Colibri observatory is currently collecting a limited amount of data while
testing the automation and processing pipeline routines. Full operational
activities are expected to begin in the summer of 2022.

## References

  * Arimatsu et al. (2017) Arimatsu, K., Tsumura, K., Ichikawa, K., Usui, F., Ootsubo, T., Kotani, T., et al. (2017).  Organized Autotelescopes for Serendipitous Event Survey (OASES): Design and performance.  _Publications of the Astronomical Society of Japan_ 69, 68.  10.1093/pasj/psx048
  * Arimatsu et al. (2019) Arimatsu, K., Tsumura, K., Usui, F., Shinnaka, Y., Ichikawa, K., Ootsubo, T., et al. (2019).  A kilometre-sized Kuiper belt object discovered by stellar occultation using amateur telescopes.  _Nature Astronomy_ 3, 301–306.  10.1038/s41550-018-0685-8
  * Bailey (1976) Bailey, M. E. (1976).  Can ‘invisible’ bodies be observed in the Solar System?  _Nature_ 259, 290–291.  10.1038/259290a0
  * Barbary (2016) Barbary, K. (2016).  Sep: Source extractor as a library.  _Journal of Open Source Software_ 1, 58.  10.21105/joss.00058
  * Bertin and Arnouts (1996) Bertin, E. and Arnouts, S. (1996).  SExtractor: Software for source extraction.  _Astronomy and Astrophysics, Supplement_ 117, 393–404.  10.1051/aas:1996164
  * Bickerton et al. (2008) Bickerton, S. J., Kavelaars, J. J., and Welch, D. L. (2008).  A Search For Sub-km Kuiper Belt Objects with the Method of Serendipitous Stellar Occultations.  _The Astronomical Journal_ 135, 1039–1049.  10.1088/0004-6256/135/3/1039
  * Bickerton et al. (2009) Bickerton, S. J., Welch, D. L., and Kavelaars, J. J. (2009).  Kuiper Belt Object Occultations: Expected Rates, False Positives, and Survey Design.  _Astronomical Journal_ 137, 4270–4281.  10.1088/0004-6256/137/5/4270
  * Chang et al. (2006) Chang, H.-K., King, S.-K., Liang, J.-S., Wu, P.-S., Lin, L. C.-C., and Chiu, J.-L. (2006).  Occultation of X-rays from Scorpius X-1 by small trans-neptunian objects.  _Nature_ 442, 660–663.  10.1038/nature04941
  * Delsanti and Jewitt (2006) Delsanti, A. and Jewitt, D. (2006).  The Solar System Beyond The Planets.  In _Solar System Update_ , eds. P. Blondel and J. W. Mason. 267\.  10.1007/3-540-37683-6_11
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., Prusti, T., de Bruijne, J. H. J., Babusiaux, C., et al. (2021).  Gaia Early Data Release 3. Summary of the contents and survey properties.  _Astronomy and Astrophysics_ 649, A1.  10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., Brown, A. G. A., Vallenari, A., Babusiaux, C., et al. (2016).  The Gaia mission.  _Astronomy and Astrophysics_ 595, A1.  10.1051/0004-6361/201629272
  * Harris (1990) Harris, W. E. (1990).  A comment on image detection and the definition of limiting magnitude.  _Publications of the Astronomical Society of the Pacific_ 102, 949–953.  10.1086/132720
  * Jones et al. (2006) Jones, T. A., Levine, A. M., Morgan, E. H., and Rappaport, S. (2006).  Millisecond Dips in Sco X-1 are Likely the Result of High-Energy Particle Events.  _The Astronomer’s Telegram_ 949, 1
  * Lang et al. (2010) Lang, D., Hogg, D. W., Mierle, K., Blanton, M., and Roweis, S. (2010).  Astrometry.net: Blind Astrometric Calibration of Arbitrary Astronomical Images.  _The Astronomical Journal_ 139, 1782–1800.  10.1088/0004-6256/139/5/1782
  * Lehner et al. (2012) Lehner, M. J., Wang, S.-Y., Alcock, C. A., Cook, K. H., Furesz, G., Geary, J. C., et al. (2012).  The Transneptunian Automated Occultation Survey (TAOS II).  In _Ground-based and Airborne Telescopes IV_ , eds. L. M. Stepp, R. Gilmozzi, and H. J. Hall. vol. 8444 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 84440D.  10.1117/12.925541
  * Lehner et al. (2009) Lehner, M. J., Wen, C. Y., Wang, J. H., Marshall, S. L., Schwamb, M. E., Zhang, Z. W., et al. (2009).  The Taiwanese-American Occultation Survey: The Multi-Telescope Robotic Observatory.  _Publications of the Astronomical Society of the Pacific_ 121, 138\.  10.1086/597516
  * Levison and Stern (2001) Levison, H. F. and Stern, S. A. (2001).  On the size dependence of the inclination distribution of the main kuiper belt.  _The Astronomical Journal_ 121, 1730–1735.  10.1086/319420
  * Pass et al. (2018) Pass, E., Metchev, S., Brown, P., and Beauchemin, S. (2018).  Pipeline for the Detection of Serendipitous Stellar Occultations by Kuiper Belt Objects with the Colibri Fast-photometry Array.  _Publications of the Astronomical Society of the Pacific_ 130, 014502\.  10.1088/1538-3873/aa971f
  * Roques et al. (2006) Roques, F., Doressoundiram, A., Dhillon, V., Marsh, T., Bickerton, S., Kavelaars, J. J., et al. (2006).  Exploration of the Kuiper Belt by High-Precision Photometric Stellar Occultations: First Results.  _The Astronomical Journal_ 132, 819–822.  10.1086/505623
  * Roques and Moncuquet (2000) Roques, F. and Moncuquet, M. (2000).  A detection method for small kuiper belt objects: The search for stellar occultations.  _Icarus_ 147, 530 – 544.  https://doi.org/10.1006/icar.2000.6452
  * Roques et al. (1987) Roques, F., Moncuquet, M., and Sicardy, B. (1987).  Stellar occultations by small bodies - Diffraction effects.  _The Astronomical Journal_ 93, 1549–1558.  10.1086/114438
  * Schlichting et al. (2012) Schlichting, H. E., Ofek, E. O., Sari, R., Nelan, E. P., Gal-Yam, A., Wenz, M., et al. (2012).  Measuring the Abundance of Sub-kilometer-sized Kuiper Belt Objects Using Stellar Occultations.  _The Astrophysical Journal_ 761, 150.  10.1088/0004-637X/761/2/150
  * Schlichting et al. (2009) Schlichting, H. E., Ofek, E. O., Wenz, M., Sari, R., Gal-Yam, A., Livio, M., et al. (2009).  A single sub-kilometre Kuiper belt object from a stellar occultation in archival data.  _Nature_ 462, 895–897.  10.1038/nature08608
  * Vida et al. (2021) Vida, D., Šegon, D., Gural, P. S., Brown, P. G., McIntyre, M. J. M., Dijkema, T. J., et al. (2021).  The Global Meteor Network - Methodology and first results.  _Monthly Notices of the Royal Astronomical Society_ 506, 5046–5074.  10.1093/mnras/stab2008
  * Zhang et al. (2013) Zhang, Z. W., Lehner, M. J., Wang, J. H., Wen, C. Y., Wang, S. Y., King, S. K., et al. (2013).  The TAOS Project: Results from Seven Years of Survey Data.  _The Astronomical Journal_ 146, 14.  10.1088/0004-6256/146/1/14


# OpenTPS - Open-source treatment planning system for research in proton
therapy

Sophie Wuyckens  Université catholique de Louvain, Institut de Recherche
Expérimentale et Clinique (IREC), Molecular Imaging, Radiotherapy and
Oncology, 1200 Woluwe-Saint-Lambert, Belgium  Damien Dasnoy  Université
catholique de Louvain, Institute of Information and Communication Technologies
(ICTEAM), Louvain-La-Neuve, 1348, Belgium  Guillaume Janssens  Université
catholique de Louvain, Institute of Information and Communication Technologies
(ICTEAM), Louvain-La-Neuve, 1348, Belgium  Ion Beam Applications SA, Louvain-
La-Neuve 1348, Belgium  Valentin Hamaide  Université catholique de Louvain,
Institute of Information and Communication Technologies (ICTEAM), Louvain-La-
Neuve, 1348, Belgium  Ion Beam Applications SA, Louvain-La-Neuve 1348, Belgium
Margerie Huet  Université catholique de Louvain, Institut de Recherche
Expérimentale et Clinique (IREC), Molecular Imaging, Radiotherapy and
Oncology, 1200 Woluwe-Saint-Lambert, Belgium  Estelle Loÿen  Université
catholique de Louvain, Institute of Information and Communication Technologies
(ICTEAM), Louvain-La-Neuve, 1348, Belgium  Gauthier Rotsart de Hertaing
Université catholique de Louvain, Institute of Information and Communication
Technologies (ICTEAM), Louvain-La-Neuve, 1348, Belgium  Benoit Macq
Université catholique de Louvain, Institute of Information and Communication
Technologies (ICTEAM), Louvain-La-Neuve, 1348, Belgium  Edmond Sterpin
Université catholique de Louvain, Institut de Recherche Expérimentale et
Clinique (IREC), Molecular Imaging, Radiotherapy and Oncology, 1200 Woluwe-
Saint-Lambert, Belgium  KU Leuven, Department of Oncology, Laboratory of
experimental radiotherapy, Leuven, 3000, Belgium  Particle Therapy
Interuniversity Center Leuven - PARTICLE, Leuven, 3000, Belgium  John Lee
Université catholique de Louvain, Institut de Recherche Expérimentale et
Clinique (IREC), Molecular Imaging, Radiotherapy and Oncology, 1200 Woluwe-
Saint-Lambert, Belgium  Kevin Souris  Université catholique de Louvain,
Institut de Recherche Expérimentale et Clinique (IREC), Molecular Imaging,
Radiotherapy and Oncology, 1200 Woluwe-Saint-Lambert, Belgium  Ion Beam
Applications SA, Louvain-La-Neuve 1348, Belgium  Sylvain Deffet  Université
catholique de Louvain, Institut de Recherche Expérimentale et Clinique (IREC),
Molecular Imaging, Radiotherapy and Oncology, 1200 Woluwe-Saint-Lambert,
Belgium

###### Abstract

Introduction. Treatment planning systems (TPS) are an essential component for
simulating and optimizing a radiation therapy treatment before administering
it to the patient. It ensures that the tumor is well covered and the dose to
the healthy tissues is minimized. However, the TPS provided by commercial
companies often come with a large panel of tools, each implemented in the form
of a black-box making it difficult for researchers to use them for
implementing and testing new ideas. To address this issue, we have developed
an open-source TPS.
Approach. We have developed an open-source software platform, OpenTPS
(opentps.org), to generate treatment plans for external beam radiation
therapy, and in particular for proton therapy. It is designed to be a flexible
and user-friendly platform (coded with the freely usable Python language) that
can be used by medical physicists, radiation oncologists, and other members of
the radiation therapy community to create customized treatment plans for
educational and research purposes.
Result. OpenTPS includes a range of tools and features that can be used to
analyze patient anatomy, simulate the delivery of the radiation beam, and
optimize the treatment plan to achieve the desired dose distribution. It can
be used to create treatment plans for a variety of cancer types and was
designed to be extended to other treatment modalities.
Significance. A new open-source treatment planning system has been built for
research in proton therapy. Its flexibility allows an easy integration of new
techniques and customization of treatment plans. It is freely available for
use and is regularly updated and supported by a community of users and
developers who contribute to the ongoing development and improvement of the
software.

_Keywords_ Proton Therapy $\cdot$ Treatment Planning $\cdot$ TPS $\cdot$
Software

##  1 Introduction

Thanks to the ballistic properties of heavy charged particles, proton therapy
has the potential to better confine the dose to the target volume and spare
healthy tissues more than conventional radiation therapy. In practice,
however, the full potential of proton therapy is limited by its maturity gap
with conventional radiation modalities. For example, the first 3D intensity
modulated proton therapy (IMPT) treatment was delivered in 1999, at a time
when volumetric arc therapy (VMAT) treatments were becoming popular in photon
therapy [1]. Similarly, while VMAT was first used on a patient in 1994 [2],
the first prototype of intensity modulated proton arc therapy was only tested
in 2019 [3]. Despite challenges, proton therapy is widely acknowledged as a
beneficial option for some patients undergoing radiation treatment. However,
the full potential of proton therapy systems cannot yet be fully realized due
to limitations in the software, particularly in the treatment planning aspect.

Over the past three decades, proton therapy has benefited from a growing body
of research aimed at realizing its full potential. This has led to the
development of many computational tools, such as REGGUI [4], MCsquare [5],
matRad, and MIROpt [6], among others. These tools are widely available and
have helped accelerate research in proton therapy. In the early 2020s, we
embarked on a project to develop our own common research platform for the
various radiation therapy projects in tha lab. Using the widely-available and
user-friendly Python programming language, the platform was designed to serve
as a flexible framework for both personal and collective research efforts.
Additionally, the platform was developed to be compatible with the increasing
significance of artificial intelligence and to facilitate the distribution of
computational tasks across multiple computers.

We propose a treatment planning system (TPS) that is specifically designed for
research, called openTPS. Its open-source code is publicly available under the
Apache-2.0 license. The code is primarily written in python and designed to be
easily extended and customized. The software is organized into two packages:
the Core package, which is a library that defines data classes, data
processing methods, and IO methods, and the GUI package, which offers a
graphical interface for viewing and interacting with the data.

The Core package of OpenTPS is the main software library and includes a range
of features that are essential for proton therapy treatment planning. Some of
the key features available in the Core package of OpenTPS include:

  1. A.

Data management and processing: importing, exporting and processing patient
data.

  2. B.

Dose computation: computing the dose from proton therapy treatment plans using
fast Monte Carlo simulations with MCsquare.

  3. C.

Treatment planning: creating treatment plans for individual patients,
including the ability to define beam angles, energy levels, spot spacing, etc.

  4. D.

Treatment evaluation: performing quality assurance checks on treatment plans,
including the ability to evaluate robustness, compare doses, and assess the
conformity of the treatment plan to the patient’s anatomy.

  5. E.

Data augmentation: creating artificial variations of patient data for
robustness evaluation and for machine learning applications.

In the next section, we detail these tools and give an example of their use
for the development of new planning methods though an introduction to the 4D
image support suite which is included in OpenTPS. Assurance tools for dose
quality are discussed in section 2.

##  2 GUI package

OpenTPS comes with an optional graphical user interface (GUI) for
visualization of 3D images, dose maps, contours and dynamic images. Several
other interactive tools are available: profile plots, visualization of DVH,
segmentation, etc. The implementation is based on the PyQt, VTK, and PyQtGraph
libraries, which are widely recognized for their capabilities in data
visualization and analysis in scientific and medical applications. Figure 1
shows a screenshot of the GUI loaded with a head-and-neck patient geometry and
the optimized treatment dose associated to the case. In addition, graphical
modules have been implemented to design and optimize a treatment plan. The
entire GUI is extensible and modifiable by the user. Examples of extensions
can be found in the OpenTPS documentation.

Figure 1: Main window of the OpenTPS GUI.

##  3 Core package

###  3.1 Data management and processing

The current version of OpenTPS supports the import patient data in DICOM
format (CT, RT structures, RT dose, RT ion plan and vector field), MHD format
(3D images), and OpenTPS-specific serialized object format. DICOM export is
currently supported for treatment plans and doses. Any images can also be
exported in MHD format.

The core library of OpenTPS includes a range of common image processing
methods: resampling, coordinate system conversion, affine transformations, and
filtering. Additionally, it has segmentation options such as automatic
segmentation for body, lung and bones in CT images.

Rigid and deformable image registration methods are available, including the
diffeomorphic Morphons algorithm [4]. This algorithm computes 3D to 3D
deformations that are consistent with the anatomy of the organs. The
deformations are represented by a velocity field and/or the more usual
displacement field. While algebraic operations - such as summing and scaling
fields - can be performed on velocity fields, they cannot be applied directly
on displacement fields. On the other hand, when a deformation is applied on 3D
data, the displacement field must be used. The displacement is automatically
computed from the velocity field when the deformation is applied using a field
exponentiation operation.

The morphons registration is also used to create motion models. Deformation
fields between the 4DCT phases are computed. Then with the deformations from
all 4DCT phases, the midposition (MidP) image [7] can be computed, as
illustrated in Figure 2. A motion model, called Dynamic3DModels in OpenTPS, is
composed of the MidP image and the various deformation fields from the MidP
image to each of the 4DCT phases.

Figure 2: The tumor motion is schematically represented by the hysteresis
formed by the tumor position (gray circles) at each phase of the 4DCT (a).
Deformation fields (blue vectors) are generated using registration between the
first phase (P1) and all others phases (b). Then, all phases can be deformed
to the average motion called the Mid- Position where the MidP-CT image is then
created by taking the median in each voxel of all deformed phases (c).

Digitaly reconstructed radiographs (DRR) can also be computed using a 3DCT to
simulate in room x-ray devices.

###  3.2 Dose computation

OpenTPS was designed with the flexibility to support multiple dose engines as
a way to provide users with a range of options for proton therapy treatment
planning. It also allows users to take advantage of the strengths of different
dose engines, such as the accuracy and precision of Monte Carlo-based engines,
or the speed and efficiency of other types of engines.

A stable release of the fast Monte Carlo named MCsquare is included in OpenTPS
[5]. It is available via the class named MCsquareDoseCalculator and can be
used to compute dose distributions, dose influence matrices, and LET
distributions. Input paramaters such as the number of primaries, the beam
model and the CT calibration are required to compute a dose based on CT and a
treatment plan. From those inputs, MCsquareDoseCalculator generates the files
required by MCsquare. Then, MCsquare simulates the interaction of each
particle in batches until it reaches a user-defined minimal statistical
uncertainty or it simulates all the primaries. Mcsquare generates the output
dose in mhd format in the OpenTPS workspace, a folder located by default in
the home directory. Figure 3 summarizes the dose calculation workflow in
OpenTPS.

Figure 3: Example of dose workflow in OpenTPS

###  3.3 Treatment planning

The treatment planning module of OpenTPS provides a set of tools and functions
that are used to create and optimize a treatment plan based on several inputs.
In general, the method followed depends on the treatment delivery modality
chosen. As OpenTPS is initially built for research in proton therapy, the
implementation of the treatment planning module conform to IMPT delivery. In
pencil beam scanning (PBS), the dose is painted using scanning magnets that
steer each pencil beam to a specific location in the target volume, namely a
spot. Dividing the tumor into iso-energy layers, the treatment is therefore
delivered by irradiating the volume spot-by-spot and layer-by-layer. This
process is repeated for each irradiation direction. The intensity of each
pencil beam being modulated, this treatment mode results in a very high dose
conformity and a large flexibility in the planning allowing for better organ
at risk sparing compared to others such as scattering techniques [8]. The
modulation of the spot intensity is the aim of plan optimization.

####  3.3.1 Plan optimization

The treatment planning module operates in several steps. The first stage
consists in creating a PlanDesign object. Loading a CTImage and RTStruct
associated to it, the planner defines the gantry and couch angles for the
given disease site. Based on this information, spots are placed in beam-eye-
view on a hexagonal grid through a ray-tracing procedure that computes min-max
water-equivalent path lengths (translated into energies) covering the target
volume and then spaced by pre-defined spot and layer spacings for each beam
direction. Internally, a RTPlan object is initialized with the angles, energy
layers and spot positions. The spot intensities are left behind for the moment
and initialized to one. Before proceeding to the optimization, the dose
distribution of each spot onto the patient voxelized-geometry, namely a
beamlet, must be computed. By default, the Monte Carlo dose engine named
MCsquare is used to simulate, for each spot, the physics of proton
interactions with the matter and score the dose deposited in each voxel.
Nevertheless, the framework is flexible enough to implement any other dose
calculation engine. The output of this step is a beamlet matrix, named
$\mathcal{A}$, which is stored in memory. To conserve memory space and enhance
computational efficiency during optimization, a sparse matrix format is
utilized. The number of rows in $\mathcal{A}$ corresponds to the number of
voxels in the patient (resampled or not) and the number of columns corresponds
to the number of spots in the plan. Finally, the clinical objective are added
to the PlanDesign object. The region of interest (ROI), the metric (min, max
or mean), the dose limit value and the objective weight must be specified for
each FidObjective. An optional parameter can be set to make the objective
robust against uncertainties.

Figure 4: Objectives panel of the OpenTPS GUI

When the plan and the associated beamlets are created, the second stage of
treatment planning can now take place. Knowing that the total dose delivered
to the patient can be computed as the sum of all beamlets weighted by the
intensity of the corresponding spots, the remaining task is to optimize the
weights $\bm{x}$ so that the total dose $\bm{d}$ meets the optimization
objectives. The optimization problem is stated as follows:

$\displaystyle\begin{split}\min_{x}f(\bm{d})\\\ \text{s.t.}\hskip
20.00003pt&\bm{d}=\mathcal{A}\bm{x}\\\ &\bm{x}\geq\bm{0}\end{split}$ (1)

The total objective function $f(\bm{d})$ is defined as a weighted sum of
asymmetric quadratic objectives with squared ramps used as surrogate to
constraint on a minimum, maximum or on a mean dose. The total objective
function has therefore usually four terms corresponding to the targets ($T$)
and the OARs ($O$), which is written as follows:

$f(\bm{d})=\sum_{t\in T}\left(\frac{1}{|v_{t}|}\sum_{i\in
v_{t}}\left(w^{+}_{t}\cdot(d_{i}-p_{t}^{max})_{+}^{2}+w^{-}_{t}\cdot(d_{i}-p_{t}^{min})_{-}^{2}\right)\right)\\\
+\sum_{o\in O}\left(\frac{w^{+}_{o}}{|v_{o}|}\sum_{i\in
v_{o}}(d_{i}-p_{o}^{max})_{+}^{2}\right)+\sum_{o\in
O}\left(w_{o}(\overline{\bm{d}}-p_{o}^{maxmean})_{+}^{2}\right)$ (2)

where $p$ is the prescribed dose and the notations $(\cdot)_{+}$ and
$\overline{\bm{\cdot}}$ represent $max\\{\cdot,0\\}$ and
$\frac{1}{|v|}\sum_{i\in v}d_{i}$ respectively.

OpenTPS comes with several built-in optimizers that solve the above problem.
It is also designed so that other optimizers may be easily implemented. A
variety of iterative solvers are provided: quasi-Newton’s methods from Scipy,
in-house implementation of LBFGS algorithms [9], gradient-based methods such
as the classical gradient descent and the FISTA algorithm [10] from
PyUNLockBox [11]. FISTA solves non-differentiable convex optimization problems
by proximal splitting. For example, it can be of particular interest for ArcPT
treatment plan optimization where one wants to maximize the spot sparsity by
adding a $l_{1}$-norm penalty to the total objective function.

Linear programming [12] methods are also available and can be used to solve
the optimization problem. Based on the list of objectives in the plan design,
a model is built with linear equations instead of the quadratic version
defined in Eq. 2. The model and the solver are built on the Gurobipy package
[13], a Python interface to the powerful Gurobi solver, which requires a
specific license. Gurobi includes a variety of algorithms to solve linear
programming problems, including the simplex algorithm [14], the interior point
method [15], and the concurrent optimizer.

Continual research is being conducted to investigate new methods, such as the
beamlet-free algorithm. This innovative approach simulates the dose in proton
batches of randomly selected spots and evaluates their impact on the objective
function in each iteration. Using an estimated gradient, the spot weights are
then updated, leading to a new spot probability distribution from which the
next spot will be sampled. This approach holds potential for further improving
the treatment planning process in terms of time and memory requirements.

The convergence rate of each of these algorithms may depend on the specific
characteristics of the optimization problem. A summary of available solvers
and their capabilities is provided in Table 1. The “order” refers to the
number of derivatives used in the optimization algorithm. First-order methods
only use the first derivative (i.e., the gradient), while second-order methods
use the first and second derivatives (i.e., the Hessian matrix).

Solver | Order | Solving function | Convergence
---|---|---|---
Type | Rate
Gradient descent | 1 | Smooth convex | Linear |  | O(1/k)
---
where $k$: # of iterations
FISTA | 1 |  | Smooth convex
---
Proximable non-smooth
Linear |  | O($1/k$) for smooth convex functions
---
O($1/k^{2}$) for proximable
non-smooth functions
| Stochastic gradient descent
---
(“Beamlet-free”)
1 | Smooth convex | / |  | Depends on the problem, parameters
---
and learning rate.
BFGS* | 2 | Smooth convex | Superlinear |  O($1/k^{2}$) where k: # of iterations
L-BFGS*
| L-BFGS-B
---
(Scipy-only)
| Smooth convex
---
with bounds constraints
LP** | Simplex | / |  Linear objective functions Linear equality and inequality constraints | / |  | O($2^{n}$) - worst case
---
where n: # of variables
Interior point | Polynomial
Concurrent | Depends on the problem
Table 1: Summary of available optimizers in OpenTPS. *Scipy and in-house
versions available. **Gurobi optimizer.

When the algorithm is selected and the plan is provided, an IMPTPlanOptimizer
object is created. Additional parameters can be provided such as the desired
stopping criterion for the optimization. The method optimize will start the
optimization using the solver supplied. The result holds, among others, the
optimized spot weights and the optimized dose computed from the beamlet
matrix. Figure 5 summarizes the different steps in the optimization workflow
in OpenTPS.

Figure 5: Optimization workflow in OpenTPS.

The structure of the plan optimization module is divided in three main sub-
modules. Firstly, the optimization engines described above are located in the
solvers module. Secondly, to compose the objective function, several functions
classes such as the classical dose objectives and more specialized objectives
(norms, projections, etc.) are grouped in the objectives module. The user can
also define a custom function for optimization by implementing an evaluation
method and either a proximal operator or gradient method. Thirdly,
acceleration schemes are included in the acceleration module such as
backtracking linesearch typically used for gradient descent or Newton’s
methods, backtracking acceleration based on a quadratic approximation of the
objective, FISTA acceleration for forward-backward solvers.

The GUI package of OpenTPS provides tools to analyze the resulting dose
distribution. The dose can directly be overlaid using a color map on the CT
image and contours to examine its homogeneity and conformity. Dose-volume
histograms can also be generated to inspect quality of the plan. Dose line
profiles can be obtained by manually drawing on the image in the viewer.
Finally, DVH metrics or important index (conformity [16] and homogeneity
indices) can be extracted from the DVH. Figure 6 shows an example of
optimization result display for a base-of-skull tumor brain scripted in the
core library. It includes a plot of the convergence of the objective function
evaluation along the iterations.

Figure 6: Robust optimized result obtained with scripting within the Core
module. Left: dose distribution. Middle: DVHs. Right: Convergence plot

####  3.3.2 Robust optimization

It is well-known that the PTV margin recipe [17] used in conventional photon
treatments does not apply well in IMPT due to the presence of proton range
uncertainties [18, 19, 20]. The alternative approach is to make use of robust
optimization [21, 22]. In the robust scheme, the dose gradients per field are
optimized so that they combine adequately even in the presence of range
uncertainties. Specific combinations of uncertainties are commonly referred to
as “scenarios”. They include patient motion and changes in anatomy resulting
in variations in the range of the proton beam. The main idea behind robust
optimization is the evaluation of the objective function in a discrete set of
scenarios. Multiple robust optimization approaches have been proposed
depending on the definition of the objective function and on the scenarios
[23, 19, 24, 25]. In OpenTPS, scenario-wise worst-case robust optimization is
implemented. It minimizes the objective function in the worst scenario, within
a certain confidence interval. The robust optimization problem, solved in
OpenTPS, can now be summarized as,

$\displaystyle\begin{split}\min_{x}\max_{s\in\mathcal{S}}f(\bm{d},s)\\\
\text{s.t.}\hskip 20.00003pt&\bm{d}=\mathcal{A}\bm{x}\\\
&\bm{x}\geq\bm{0}\end{split}$ (3)

where the max-operator applies to the scenarios $\\{s\in\mathcal{S}\\}$.

The current implementation allows the selection of scenarios in the error
space. The planner specifies the setup error (systematic and random) in
millimeter as well as the range uncertainties in percentage. A total of 21
scenarios will be simulated (7 setup x 3 range). In terms of workflow, after
the robustness settings are entered, the TPS will simulate the beamlets for
each scenario. Next, the planner must specify among the objective list which
ones are robust for the optimization. With that information, the optimizer
will solve the problem by updating the weights using the objective function
gradient calculated from the worst-case scenario at each iteration. The robust
optimized treatment plan produces the best possible outcome according to the
worst case scenario and robustness parameter settings.

###  3.4 Treatment evaluation

Robust evaluation involves the examination of an optimized treatment plan
under a variety of different scenarios that could occur during treatment
delivery. This approach may be used to ensure that the treatment plan is
effective and reliable under a range of possible conditions.

####  3.4.1 Robustness to setup and range uncertainties

In OpenTPS, a dose-volume evaluation is performed to assess the quality of the
treatment plan. Two strategies are available depending on the scenarios
definition. The first method, which is based on good practice rules and is
commonly used in clinical routine, involves verifying the dose in the error
space for a small number of extreme scenarios [26]. In OpenTPS, the dose is
re-computed with MCsquare for 81 scenarios resulting from a combination of the
setup and range scenarios previously described. The second strategy generates
a confidence interval in the so-called dosimetric space. It refers to all
possible realizations of the treatment plan, or the different dose
distributions that could be achieved. Described in [27], the scenarios are
randomly sampled by MCsquare. The analysis in the dosimetric space reports the
percentage of evaluation scenarios that meet a certain criterion such as
minimum dose coverage for a target percentage (e.g., 95% coverage for 90% of
scenarios).

In OpenTPS, to proceed to robust evaluation, a Robustness object (associated
to a plan) must be created and configured for the intented goal. A MCsquare
simulation will compute the robust scenarios based on these information. The
scenarios can be analyzed following the two strategies previously described
using the member methods analyzeErrorSpace or analyzeDosimetricSpace. The DVHs
of the nominal scenario and DVH bands can be recomputed for a list of ROIs and
plotted for better visualization of the robustness of the generated treatment
plan.

####  3.4.2 Robustness to intra-fraction motion (4D simulation)

Two types of dose deterioration arise while treating moving targets in proton
therapy using pencil beam scanning:

Motion-induced range variations

The tumor motion can lead to a density variation in the beam path, thereby
shifting the expected proton range and deteriorating the overall dose
distribution.

Interplay effect

The tumor and scanning beam motions happen at the same time scale, leading to
an interplay effect that can lead to hot and cold spots in the target.

To analyze those deteriorating effects, one can simulate the treatment
delivery via Monte Carlo (MC) simulations on the 4DCT. In OpenTPS, those MC
simulations are performed with MCsquare.

##### 4D dose simulation

To study the impact of the proton range variation, one can perform a 4D dose
(4DD) simulation, which consists of independently simulating the dose delivery
on each motion phase of the 4DCT and computing the average of those doses.
Thus, each motion phase is considered static and contributes equally to the
total dose. Mathematically the 4D dose is given by

$\displaystyle\mathbf{F}_{i}$
$\displaystyle=\mathrm{DIR}(\mathbf{CT}_{i},\mathbf{CT}_{\text{MidP}})\text{
for all $i=1,...,P$}$ (4) $\displaystyle\mathbf{d}_{tot}$
$\displaystyle=\frac{1}{P}\sum_{i=1}^{P}\mathbf{d}^{\text{stat}}_{i}(\mathbf{v}+\mathbf{F}_{i}(\mathbf{v}))\text{
for all $\mathbf{v}\in\mathbb{R}^{3}$}$ (5)

where $v$ is a voxel of the image, $P$ is the number of motion phases,
$\mathbf{d}^{\text{stat}}_{i}$ is the dose resulting from the treatment plan
simulation on the 3DCT of motion phase $i$, $\mathrm{DIR}(\cdot)$ is the
deformable image registration operator, $\mathbf{F}_{i}$ is the deformation
field resulting from the registration, $\mathbf{CT}_{i}$ is planning CT image
$i$ and $\mathbf{CT}_{\text{MidP}}$ is the Mid-position CT. This resulting
dose allows to analyze the impact of the proton range variation when motion is
present and isolating it from the interplay effect.

##### 4D dynamic dose simulation

To study the impact of the interplay effect, one can perform a 4D dynamic dose
(4DDD) simulation, which consists of simulating the delivery of beam spots one
by one, via a time delivery model, on a sequence of 3DCTs. The sequence of
3DCTs is generally a looped version of the planning 4DCT, but it could be a
richer set of synthetic CTs representing multiple breathing motions. The time
delivery model is a model that associates a timing with the delivery of each
spot in the treatment plan. This model is generally machine-specific and gives
an expected delivery timing depending on several parameters such as the delay
between the delivery of two adjacent spots, the time required for switching
energy, etc. Knowing the delivery timing of each spot in the treatment plan,
we can create a breathing model and know exactly in which motion phase each
spot will be delivered. The breathing model has two parameters in the simplest
case: the breathing period and the number of motion phases of the 4DCT.
Mathematically the total dose is given by

$\displaystyle\mathbf{F}_{i}$
$\displaystyle=\mathrm{DIR}(\mathbf{CT}_{i},\mathbf{CT}_{\text{MidP}})\text{
for all $i=1,...,P$}$ (6) $\displaystyle\mathbf{d}_{tot}$
$\displaystyle=\sum_{i=1}^{P}\mathbf{d}^{\text{dyn}}_{i}(\mathbf{v}+\mathbf{F}_{i}(\mathbf{v}))\text{
for all $\mathbf{v}\in\mathbb{R}^{3}$.}$ (7)

where $\mathbf{d}^{\text{dyn}}_{i}$ is the dose associated with the motion
phase $i$ that includes all spots predicted to be delivered in that motion
phase. Repeating the simulation by starting the delivery at different motion
phases allows for quantifying the uncertainty of the final dose distribution.
Note that we suppose each motion phase is a static image, thus neglecting
motion during a motion phase.

###  3.5 Data augmentation

Moving targets, such as lung or liver tumors, come with an additional
challenge: intra-fraction motion such as breathing, and usually higher inter-
fraction motion than static target such as baseline shifts. To simulate
treatment delivery on this kind of tumors, tools have been added to OpenTPS to
add intra- and inter-fraction motions, as well as creating training and test
data based on motion models (see section 3.1) to evaluate innovative methods
or AI models.

##### Inter-fraction motion simulation

Four types of inter-fraction motion can be artificially generated.

Baseline shifts

Artificial baseline shifts can be applied on 3DCT or Dynamic3DModels based on
the ROIMask of the target or organ to be translated. The baseline shift is
constructed in the form of a diffeomorphic displacement vector field modeling
a local shift of the ROIMask while preserving surrounding bony structures.

Organ/target shrinks

Organ or target shrinks can be applied on Image3D or Dynamic3DModel objects,
using a ROIMask for the organ that must be shrinked. First, the ROIMask is
dilated by 1 voxel and eroded by a number of pixels depending on the user
input in mm. From these two new masks, the dilated and eroded bands
corresponding to the difference with the original mask are computed. Then, one
by one, every voxel of the eroded band is replaced by a random sample picked
from a normal distribution $N(\mu,\sigma^{2})$ where $\mu$ is the average
value of the 10 voxels from the dilated band closest to the voxel to fill.

Rotations

Rotations around the 3 main axis can be applied to Image3D or Dynamic3DModel
objects. The rotation is applied around an axis at the center of the image and
not around the image origin. Note that in case of multiple rotations, the
order is important as rotations are not commutative.

Translations

Translations in the 3 main axis can be applied to Image3D or Dynamic3DModel
objects.

An example of synthetic inter-fraction modification on a 3DCT image can be
seen in Figure 7.

Figure 7: Example of inter-fractions transforms applied to a 3DCT. Small
translations and rotations were applied to the entire image, and additionally
a baseline shift and a shrink were applied to the lung tumor.

##### Intra-fraction breathing simulation

Using a motion model, a sequence of images following a specific breathing
pattern can be created. First, synthetic 1D or 3D breathing signals can be
generated with different parameters to change the breathing period, the noise
or to add irregularities such as apneas (Figure 8).

Figure 8: Example of synthetic breathing signals, regular in blue and
irregular in orange.

Then, one or multiple points can be selected in the midP image of a motion
model to apply these breathing signals in a specified direction. If multiple
points and signals are used (for example to partially decorrelate the motion
of the tumor and the motion of the skin), they are combined using weight maps
and linear interpolations [28] to create 3D image sequences containing any
type of breathing patterns. These image sequences can be used to evaluate
treatment robustness, new treatment methods [29] or create additional patient
data to train AI methods [30].

##  4 Perspectives

OpenTPS is a unique platform for research in proton therapy treatment
planning, being an open-source tool written in a user-friendly programming
language. This allows for both internal contributions and international
collaboration to enhance the platform quickly, making it increasingly
sophisticated over time. Our research group is currently focusing on several
exciting areas within OpenTPS, including arc proton therapy treatment planning
[31, 32] and FLASH proton therapy[33]. Depending on the needs of the users,
new features might be added in the future such as (but not limited to) other
treatment or image modalities, objects tracking or position prediction in
dynamic sequences. Additionally, the platform has potential for real-time
treatment delivery visualization and 2D to 3D image reconstruction. As the
field of proton therapy continues to evolve, OpenTPS offers a flexible and
dynamic platform for exploring new possibilities and pushing the boundaries of
what is possible.

The information to get started are available on the
website111http://opentps.org/ or on
gitlab222https://gitlab.com/openmcsquare/opentps.

##  5 Acknowledgments

During the development of OpenTPS, Dr. Sylvain Deffet received funding of the
Walloon Region of Belgium through technology innovation partnership no. 8341
(EPT-1 – Emerging Proton Theraexact Phase 1) co-led by MecaTech and BioWin
clusters. Sophie Wuyckens and Kevin Souris are funded by the Walloon Region as
part of the Arc Proton Therapy convention (Pôles Mecatech et Biowin). Dasnoy
Damien, Gauthier Rotsart de Hertaing and Valentin Hamaide are supported by the
Walloon Region, SPWEER Win2Wal program project 2010149. Estelle Loÿen is a
Televie grantee of the Fonds de la Recherche Scientifique - F.N.R.S.
Computational resources have been provided by the supercomputing facilities of
the Université catholique de Louvain (CISM/UCL) and the Consortium des
Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI) funded
by the F.R.S.-FNRS under convention 2.5020.11.

## References

  * [1] K. Hatano, N. Tohyama, T. Kodama, N. Okabe, M. Sakai, and K. Konoeda,  Current status of intensity-modulated radiation therapy for prostate cancer: History, clinical results and future directions,  International Journal of Urology 26, 775–784 (2019).
  * [2] T. R. Mackie, T. Holmes, S. Swerdloff, P. Reckwerdt, J. O. Deasy, J. Yang, B. Paliwal, and T. Kinsella,  Tomotherapy: A new concept for the delivery of dynamic conformal radiotherapy,  Medical Physics 20, 1709–1719 (1993).
  * [3] X. Li, G. Liu, G. Janssens, O. De Wilde, V. Bossier, X. Lerot, A. Pouppez, D. Yan, C. Stevens, P. Kabolizadeh, and X. Ding,  The first prototype of spot-scanning proton arc treatment delivery,  Radiotherapy and Oncology 137, 130–136 (2019).
  * [4] G. Janssens, L. Jacques, J. O. D. Xivry, X. Geets, and B. MacQ,  Diffeomorphic registration of images with variable contrast enhancement,  International Journal of Biomedical Imaging (2011).
  * [5] K. Souris, J. A. Lee, and E. Sterpin,  Fast multi-purpose Monte Carlo simulation for proton therapy using multi- and many-core CPU architectures,  Medical Physics 1700, 1–23 (2016).
  * [6] A. M. Barragán Montero, K. Souris, D. Sanchez-Parcerisa, E. Sterpin, and J. A. Lee,  Performance of a hybrid Monte Carlo-Pencil Beam dose algorithm for proton therapy inverse planning,  Medical physics 45, 846—862 (2018).
  * [7] J. W. H. Wolthaus, J.-J. Sonke, M. van Herk, and E. M. F. Damen,  Reconstruction of a time-averaged midposition CT scan for radiotherapy planning of lung cancer patients using deformable registration,  Medical Physics (2008).
  * [8] H. Paganetti,  Proton therapy physics,  CRC press, 2018.
  * [9] R. Byrd, P. Lu, J. Nocedal, and C. Zhu,  A limited memory algorithm for bound constrained optimization,  SIAM Journal of Scientific Computing 16, 1190–1208 (1995).
  * [10] A. Beck and M. Teboulle,  A fast iterative shrinkage-thresholding algorithm for linear inverse problems,  SIAM Journal on Imaging Sciences 2, 183–202 (2009).
  * [11] M. Defferrard, R. Pena, and N. Perraudin,  PyUNLocBoX: Optimization by Proximal Splitting.
  * [12] G. B. Dantzig,  Linear programming,  Operations research 50, 42–47 (2002).
  * [13] Gurobi Optimization, LLC,  Gurobi Optimizer Reference Manual, 2022.
  * [14] G. B. Dantzig,  Origins of the Simplex Method, page 141–151,  Association for Computing Machinery, New York, NY, USA, 1990.
  * [15] F. A. Potra and S. J. Wright,  Interior-point methods,  Journal of Computational and Applied Mathematics 124, 281–302 (2000),  Numerical Analysis 2000. Vol. IV: Optimization and Nonlinear Equations.
  * [16] E. Shaw, R. Kline, M. Gillin, L. Souhami, A. Hirschfeld, R. Dinapoli, and L. Martin,  Radiation therapy oncology group: Radiosurgery quality assurance guidelines,  International Journal of Radiation Oncology*Biology*Physics 27, 1231–1239 (1993).
  * [17] M. van Herk, P. Remeijer, C. Rasch, and J. V. Lebesque,  The probability of correct target dosage: dose-population histograms for deriving treatment margins in radiotherapy,  International Journal of Radiation Oncology*Biology*Physics 47, 1121–1135 (2000).
  * [18] A. J. Lomax,  Intensity modulated proton therapy and its sensitivity to treatment uncertainties 1: the potential effects of calculational uncertainties,  Physics in Medicine & Biology 53, 1027 (2008).
  * [19] J. Unkelbach, T. Bortfeld, B. C. Martin, and M. Soukup,  Reducing the sensitivity of IMPT treatment plans to setup errors and range uncertainties via probabilistic treatment planning,  Medical Physics 36, 149–163 (2009).
  * [20] E. Engwall, A. Fredriksson, and L. Glimelius,  4D robust optimization including uncertainties in time structures can reduce the interplay effect in proton pencil beam scanning radiation therapy,  Medical Physics 45, 4020–4029 (2018).
  * [21] F. Albertini, E. B. Hug, and A. J. Lomax,  Is it necessary to plan with safety margins for actively scanned proton therapy?,  Physics in Medicine and Biology 56, 4399–4413 (2011).
  * [22] A. Fredriksson,  A characterization of robust radiation therapy treatment planning methods—from expected value to worst case optimization,  Medical Physics 39, 5169–5181 (2012).
  * [23] J. Unkelbach, T. C. Y. Chan, and T. Bortfeld,  Accounting for range uncertainties in the optimization of intensity modulated proton therapy,  Physics in Medicine & Biology 52, 2755 (2007).
  * [24] D. Pflugfelder, J. J. Wilkens, and U. Oelfke,  Worst case optimization: a method to account for uncertainties in the optimization of intensity modulated proton therapy,  Physics in Medicine & Biology 53, 1689 (2008).
  * [25] A. Fredriksson, A. Forsgren, and B. Hårdemark,  Minimax optimization for handling range and setup uncertainties in proton therapy,  Medical Physics 38, 1672–1684 (2011).
  * [26] E. Sterpin, S. T. Rivas, F. V. den Heuvel, B. George, J. A. Lee, and K. Souris,  Development of robustness evaluation strategies for enabling statistically consistent reporting,  Physics in Medicine & Biology 66, 045002 (2021).
  * [27] K. Souris, A. Barragan Montero, G. Janssens, D. Di Perri, E. Sterpin, and J. A. Lee,  Technical Note: Monte Carlo methods to comprehensively evaluate the robustness of 4D treatments in proton therapy,  Medical Physics 46, 4676–4684 (2019).
  * [28] D. Dasnoy-Sumell, A. Aspeel, K. Souris, and B. Macq,  Locally tuned deformation fields combination for 2D cine-MRI-based driving of 3D motion models,  Physica Medica 94, 8–16 (2022).
  * [29] V. Hamaide, K. Souris, D. Dasnoy, F. Glineur, and B. Macq,  Real-time image-guided treatment of mobile tumors in proton therapy by a library of treatment plans: a simulation study,  Medical Physics (2022).
  * [30] E. Loyen, D. Dasnoy-Sumell, and B. Macq,  3DCT reconstruction from single X-ray projection,  in 4D treatment workshop on particle therapy, 2021.
  * [31] S. Wuyckens, M. Saint-Guillain, G. Janssens, L. Zhao, X. Li, X. Ding, E. Sterpin, J. A. Lee, and K. Souris,  Treatment planning in arc proton therapy: Comparison of several optimization problem statements and their corresponding solvers,  Computers in Biology and Medicine 148, 105609 (2022).
  * [32] S. Wuyckens, L. Zhao, M. Saint-Guillain, G. Janssens, E. Sterpin, K. Souris, X. Ding, and J. A. Lee,  Bi-criteria Pareto optimization to balance irradiation time and dosimetric objectives in proton arc therapy,  Physics in Medicine & Biology 67, 245017 (2022).
  * [33] S. Deffet, K. Souris, and E. Sterpin,  Optimization of patient-specific range modulators for conformal FLASH proton therapy,  Medical Physics (submitted).