electronCT - an imaging technique using very-high energy electrons

The electronCT technique is an imaging method based on the multiple Coulomb scattering of relativistic electrons and has potential applications in medical and industrial imaging. It utilizes a pencil beam of electrons in the very high energy electron (VHEE, 50–250 MeV) range and a single detection layer for the determination of the beam profile. The technique constitutes a projectional, two-dimensional imaging method and thus also qualifies for the tomographic reconstruction of samples. Given the simplicity of the technical setup and its location behind the sample, the electronCT technique has potential synergies with VHEE radiotherapy, making use of the same electron source for both treatment and diagnostics and thus being a candidate for in situ imaging and patient localization. At the same time, several technical challenges arise from the measurement technique when applied for the imaging of living beings. Measurements performed at the ARES linear particle accelerator at an electron energy of 155 MeV using a mouse phantom and a Timepix3 silicon pixel detector assembly demonstrate the feasibility of this technique. Both projectional and tomographic reconstructions are presented and the potential and limits of the technology are discussed.

1 Introduction

In the past years, the field of radiation therapy has seen fast developments, driven by particle accelerator technologies that leverage the use of the very high energy electron (VHEE) regime, which typically refers to electrons in the energy range of 50 MeV–250 MeV, and the advancements in FLASH radiotherapy [1]; [2]; [3]; [4]. In order to achieve a safe and effective treatment in radiation therapy, reliable imaging strategies are inevitable. The application of imaging techniques in the planning and performance of radiation therapy is condensed in the term of image guided radiation therapy (IGRT). A large number of imaging methods can be and are applied in the context of IGRT, such as X-ray imaging and CT, magnetic resonance imaging (MRI), positron emission tomography (PET), ultrasound imaging and camera-based imaging [5]. While a good accuracy of a few millimeters can be achieved in aligning the treatment beam with an image using separate devices for these tasks, a unification of the treatment and an imaging device could improve this accuracy significantly and has been achieved for photon therapy [6].

We present an alternative imaging technique called electron computed tomography (electronCT), that is based on highly energetic charged particles such as VHEE and their interaction with matter, in particular the multiple Coulomb scattering of the particles traversing the patient. This measurement technique enables the use of the same accelerating structure for both treatment and diagnostics in the context of radiotherapy with VHEE and thus naturally creates a common coordinate system and does not require any alteration of the instrumentation around the patient. This makes electronCT a candidate for an IGRT imaging technique in VHEE radiotherapy, for example, for the in situ localization of tumors or as an input for the alignment of high-resolution X-ray, CT or MRI images with the coordinate system applied for the treatment.

This article presents the concepts of the electronCT technique and shows proof-of-concept measurements. Furthermore, the potential and limits of this method are discussed.

2 The electronCT technique

The technique of electronCT relies on the multiple Coulomb scattering of highly energetic particles in matter. When traversing matter, charged particles, typically of momenta in the order of few hundreds of MeV, are stochastically deflected by the electrostatic force of the material’s nuclei, leading to an effective deflection when traversing an object. The effective deflection depends on the radiation length $X_0$ and the thickness $l$ of the material, as well as on the charge and momentum $p$ of the incident particle. This scattering process and with it the effective deflection angle distribution are theoretically described by Molière’s theory [7]; [8]. The central part of this distribution of deflection angles can be described by a Gaussian distribution centered around zero. The width ${\theta }_0$ of this distribution is commonly approximated via a formula introduced by Highland [9], with the parameters revised by Lynch and Dahl [10]:

${\theta }_0=\frac{13.6\text{\hspace{0.17em}MeV}}{\beta cp}z\sqrt{\frac{l}{X_0}}\left(1+0.038\ln \left(\frac{l}{X_0}\right)\right),(1)$

with $\beta$ as the velocity in fractions of the speed of light $c$ and $z$ the charge number of the incident particle. The width of the distribution depends on the radiation length of the material traversed, increases with the thickness of the material and is reduced for higher momenta.

In the electronCT technique, the amount of material traversed by a beam is measured by determining the opening angle of a collimated beam of particles after traversing a sample. The amount of material is defined as the material’s thickness normalized to its radiation length, $ϵ=l/X_0$, and often referred to as the material budget. This measurement can be accomplished by detecting the transverse beam profile using a single detection layer downstream of the sample, employing position sensitive radiation detectors. Assuming an incident beam of charged particles with low transverse size and low divergence, the width of the transverse beam profile at a given distance behind the sample is thus a measure for the material budget along the path of the particles.

While many types of radiation detectors are applicable and should be chosen depending on the expected size and intensity of the beam delivered by the particle accelerator, this publication presents the use of silicon pixel detectors which are typically applied in high energy physics. Silicon detectors measure the amount of energy deposited in a sensitive volume via ionisation processes and a segmentation of the sensor in pixels allows for the retrieval of two-dimensional information on a particle’s traversal position [11]. They are typically optimized for the tracking of individual particles through a multi-detector setup, thus their dynamic range is optimized for a small number of particles per readout cell. In this case, silicon pixel detectors come with the benefit of less signal per readout channel and thus a lower required dynamic range compared to sensors segmented in strips.

Figure 1 shows the visualisation of the simulated acquisition of a single data point with this technique, performed via the semiconductor detector simulation framework ${\mathrm{A}\mathrm{l}\mathrm{l}\mathrm{p}\mathrm{i}\mathrm{x}}^2$ [12], which utilizes the software toolkit Geant4 [13] for the simulation of the interaction of particles with matter. The incident beam, of which the individual particle trajectories are shown in red, is scattered at the sample (blue) before it is detected in the silicon detector (grey). In addition, a small number of Bremsstrahlung photons is generated (light green). The sample in this case consists of two nested cylinders with radii in the range of small medical samples and radiation lengths in the range of tissues. Figure 2 shows two examples for simulated beam profiles, with (right) and without sample present (left). Both beam profiles are displayed as the charge collected in the sensor, in units of kiloelectrons, as a function of the impact pixel coordinates and can clearly be distinguished from each other. Gaussian distributions are fitted to the projections of the beam profiles onto the $x$- and the $y$-axes and the corresponding beam sizes are determined as widths of the fitted distributions in units of pixels. The means of these widths, ${\sigma }_b=\left({\sigma }_x+{\sigma }_y\right)/2$, are indicated in the corresponding graphs in Figure 2. The simulation demonstrates that a measurement of the beam profile of an initially collimated electron beam via a pixelated silicon sensor is sensitive to the traversal through a sample prior to the detection. In addition, this simulation was applied for designing the experimental layout. For this simulation, Allpix Squared version 3.0.3 [14], compiled with Geant4 version 11.2.1, has been used to simulate the interaction of a beam containing 1,000 primary electrons with the sample and the detector, applying a Geant4 physics list with the electromagnetic constructor Livermore. It should be noted that for the visualisation in Figure 1A beam of only one hundred electrons has been simulated for a better visibility.

Figure 1

Figure 1. Visualisation of the simulation corresponding to the acquisition of a single data point acquired via the electronCT technique. A beam consisting of one hundred electrons with a width of 100 µm (red lines) is scattered at a sample (blue) and consecutively passes through a silicon detector (grey). Along the electron trajectories, a small number of Bremsstrahlung photons is generated (light green).

Figure 2

Figure 2. Simulation of beam profiles as detected by the silicon pixel detector in a simulation such as shown in Figure 1 for the beam not passing through the sample (left) and traversing the center of the sample (right). Shown is the collected charge in units of kiloelectrons as a function of the impact pixel coordinates.

Two-dimensional or projectional imaging can be achieved by the bunched particle beam rastering the sample or by moving the sample across a fixed beam, such that each electron bunch and thus each measured width of a bunch profile can be attributed to a certain impact position on a virtual transverse plane through the sample.

Three-dimensional imaging is enabled by recording projections at different impact angles. This can be achieved by either rotating the sample or, similarly to X-ray based computed tomography (CT) measurements, the rotation of the particle source, here the accelerating structure, and detector around the sample. The rotation of the particle accelerator around the sample or patient is technologically available in a few modern radiation treatment facilities [15]; [16].

Consequently, the experimental setup comprises an accelerating structure delivering a beam with momenta of a few hundreds of MeV and a beam size in the order of a few 100 μm, and a silicon pixel detector capable of coping with a high data rate.

A similar technique, also based on multiple Coulomb scattering but applying the measurement of individual electron trajectories of GeV-electrons in a large beam, has proven to provide 3D imaging with resolutions in the order of 100 µm and good contrast-to-noise ratios for a wide range of material densities, at the downside of a low particle rate and thus extensive measurement times [17]. The development of scanning the sample with a pencil beam brings the opportunity to drastically decrease the measurement time. In addition, it has the potential to reduce the complexity of both the setup and the analysis by omitting the need to reconstruct individual particle trajectories and thus to reduce the time required for the reconstruction of an image.

3 Experimental setup

The measurements shown in this work were performed at the Accelerator Research Experiment at SINBAD (ARES) [18] at DESY, Hamburg, using a Timepix3 silicon pixel detector assembly [19]. A medical phantom resembling a mouse was placed between the electron extraction window of the accelerator and the detector. The individual components were positioned as closely together as possible while assuring a clearance of the motion stages with respect to further components, leading to distances between the beam window and the phantom of 68 mm and between the beam window and the silicon detector of 134 mm. A picture of the main experimental setup can be seen in Figure 3.

Figure 3

Figure 3. Picture of the electronCT setup at the ARES accelerator facility. The mouse phantom is mounted on a positioning stage and located in the center of the experiment. The beam exit window, secured by a black metallic cover, is to the left and the detector on the right hand side.

The measurements presented herein use a sample mounted to a 4D positioning system which rasters the sample across a fixed beam.

3.1 ARES

The ARES accelerator is a linear electron accelerator delivering ultra-short bunches with a charge of up to few hundreds of picocoulombs [18]. It is designed for a kinetic energy of up to 155 MeV [18] and a bunch repetition rate of between 1 Hz and 50 Hz. For this experiment a repetition rate of 10 Hz was applied and significantly lower bunch charges, estimated to be in the order of 1 fC and less, were achieved via the use of the dark current, a current generated by field emission in the accelerating structure [20]. The beam was transversely focussed onto the sample using the last quadrupole magnets of the accelerator beamline. With this, a symmetric beam spot with an RMS size of around 320 µm was reached at the position of the sample and of about 360 μm at the detection plane, as will be presented below. The transverse beam size was observed to be dominated by the scattering at the beam exit window, consisting of a titanium (Ti Grade 5) foil of 50 (5) µm thickness.

3.2 Timepix3 silicon pixel detector

The choice of the detector is a crucial parameter for the electronCT method. For the presented application, a large readout buffer and a large dynamic range are two main demands on the detection layer. The former is motivated by the fact, that at the traversal of a single bunch a large area of the detector and thus a large quantity of pixels is hit within a small amount of time, and all buffered data has to be read out before the arrival of the succeeding bunch to avoid potential dead times. The requirement for a large dynamic range arises from the fact, that the number of incident electrons remains constant, while the width of the beam profile strongly changes depending on the traversed material. In this experiment, both large peaks and long tails in the beam profile need to remain resolvable.

One type of detectors that fulfil these requirements are assemblies based on the Timepix3 or Medipix3 readout chips, which are already used in several medical applications [21]; [22]. Timepix3 readout chips feature an array of 256 × 256 pixels with a pitch of 55 µm $\times$ 55 µm and provide charge information via Time-over-threshold (ToT) detection with a TDC resolution of 10 bit, as well as a time-of-arrival (ToA) measurement. The assembly used herein consists of a Timepix3 readout chip bump-bonded to a silicon sensor with a thickness of 100 µm. The sensor was operated at a bias voltage of $-21\mathrm{V}$ to ensure full depletion [23] and read out using a Katherine readout system [24] controlled via the TrackLab software framework [25].

3.3 Phantom and positioning

The phantom represents a mouse body with the tissue additively manufactured from resin, a skeleton made from a water-gypsum mixture and organs formed from agarose [26]. These materials have been chosen due to their manufacturing properties, low costs and their X-ray attenuation coefficients close to those of a real mouse. It should be noted that the radiation length, as the underlying material parameter for electronCT measurements, does not scale linearly with the X-ray attenuation coefficient, but both strongly depend on the atomic number. However, the X-ray attenuation coefficient also exhibits a dependency on the particle energy, which is not the case for the radiation length. Hence, although the materials contained in the sample are not expected to reproduce the exact multiple Coulomb scattering behaviour of a biological sample, it is expected to provide a good benchmark for electronCT measurements.

Linear translation stages of the type LTS300 by Thorlabs [27]. were utilized for scanning the phantom across the beam in the two transverse dimensions, while the beam was kept at a constant position and the detector mounted on a fixed stand. The stages feature an accuracy of about 2 µm. The sample can also be moved along the beam axis, but the range of motion in this dimension is severely restricted as the experiment setup is optimized for minimum distances along the beam axis. The sample is furthermore mounted to a rotational stage of the type PRM1/MZ8 by Thorlabs [28]. with a vertical rotation axis, featuring a sub-degree rotational precision. This allows for an illumination from different angles and thus a tomographic measurement.

4 Methodology

4.1 Scan sequence

For two-dimensional measurements, the sample is moved across the $x$-$y$-plane, with the $x$-axis as the horizontal, and the $y$-axis as the vertical axis. The scans have been conducted using a continuous motion along the $x$-axis at discrete steps along the $y$-axis, forming a serpentine path. In this schematic, the minimum achievable resolution of the resulting images is given by the velocity of the motion stage along $x$ and the bunch repetition rate of 10 Hz, and the step size along $y$. With a motion velocity of e.g., 5 mm s⁻¹, the bunches sample the phantom with a spacing of $\mathrm{\Delta }x=$ 0.5 mm.

Three-dimensional, tomographic measurements via electronCT require a rotation of the sample, which is achieved by a rotation around the $y$-axis with the angle denoted as $\phi$. For these studies, the motion system was set up to perform two-dimensional scans as mentioned above for a configurable number of rotation angles consecutively.

An alternative sequence has been defined as a scan in the $x$-$\phi$-space at a fixed $y$-position, representing a single horizontal line for several rotation angles, which allows for the tomographic reconstruction of a single slice of the sample through the $x$-$z$-plane. The advantage of this measurement in comparison with a full three-dimensional scan is the potential to study the tomographic reconstruction potential and performance at a drastically reduced measurement time.

It should be noted that for all scan types the sequence of data points taken is not relevant for the data analysis and can be optimized e.g., for reducing the measurement time.

4.2 Data acquisition

For electronCT measurements, the electron bunch is focused at the detector, such that in most cases several electrons contribute to the signals of individual pixels. The time structure of the bunch does not allow for a separate detection of the individual primary particles and hence the summed deposited charge below each pixel is measured. As a consequence, a high bunch charge focused onto a small number of pixels would lead to a saturation of the detector front-end dynamic range and should thus be avoided.

The Timepix3 chip was configured with a threshold of about four kiloelectrons. It should be noted that the unit kiloelectrons (ke⁻) here denotes the number of charge carriers collected per pixel and not the number of primary electrons contributing to the signal. The threshold corresponds to about 40% of the signal induced by a single electron and is thus sufficient for the detection of individual beam particles.

The readout chip performs a zero-suppression, i.e., only pixels detecting a signal larger than the configured threshold are registered, and can be operated in two modes, the so-called sequential (often called frame-based) readout mode or the data-driven mode. In the frame-based mode a readout frame is defined by an external signal and the data from all pixels that registered a signal above the configurable threshold within this frame is read out after the frame has ended. The Katherine readout system allows to delay the frame start with respect to an external signal and to set a configurable frame duration. In this readout mode the system can be configured, such that each frame represents the signals from an individual bunch. The data-driven readout mode represents a continuous detector readout, leading to a data set of all signals exceeding the threshold, containing the corresponding pixel addresses including timestamps. In this mode, the timestamps of the signals can be used to group them into signals arising from individual bunches during the post-processing of the data.

In case of the frame-based readout, frames have been triggered via the accelerator machine clock, with the delay manually adjusted, such that bunches arrive within the first microsecond of a frame. The frames were configured to feature a width of 10 µs The data contain the coordinates and ToT values of all pixels registering a signal above the threshold per readout frame as well as their ToA within the frame.

In the data-driven readout mode, the data contain the coordinates and ToT values of all signals above the threshold as well as each signal’s ToA within the data stream. While the global timestamp of the readout chip itself features a range of 409.6 µs, the Katherine DAQ system is capable of extending the timestamp to a range of more than 20 days and is hence sufficient for all measurements performed.

4.3 Data analysis

4.3.1 Data processing

The data were converted and every frame was interpreted using the software framework Corryvreckan [29]. To interpret data recorded with the Katherine readout system via this framework, a module for reading data recorded by TrackLab was added. In the case of measurements applying the data-driven readout mode, no inherent frame structure is available that would allow to assign the detector data to a given state of the motion system. Instead, a continuous array containing all pixels with a signal above the threshold with the corresponding ToA is obtained. An example of the time structure of the recorded data is shown in Figure 4 by means of the number of hits as a function of their corresponding time stamps. The time structure shows a clear, regular pattern caused by the ARES bunch structure leading to large numbers of hit pixels in a vanishingly small time frame with a distance of 100 ms as expected from the bunch repetition rate of 10 Hz. Infrequent signals are recorded in between the arrivals of two bunches as seen at the positions of 12.8 s and 13.3 s and can be caused by noise or uncorrelated radiation.

Figure 4

Figure 4. Time structure of the beam generated by ARES as recorded by the Timepix3 detector assembly, represented as the number of hit pixels as a function of time.

The ToA information of each signal was used to group them into frames, where every frame represents the data induced by a single bunch. This was achieved with a minimal bias on the actual event separation by an algorithm splitting the data stream at positions where no significant group of signals has been recorded.

Figure 5 shows two charge maps for individual frames, representing a mapping of all hits to their pixel coordinates with the ToT as color scale, one corresponding to a frame in which the beam does not pass through the sample (left) and one to a frame for which the sample was traversed by the beam (right). It is readily visible that the multiple scattering of the beam particles at the object enlarges the beam laterally at the position of the sensor. The sensor front-end saturates in case of large amounts of particles impinging in individual pixels, leading to a plateau visible in Figure 5 (left).

Figure 5

Figure 5. Charge map of the detector signal from two individual bunches, one not traversing the sample (left) and one traversing it (right).

For the analysis of this data, a dedicated module in the Corryvreckan data analysis framework projects these maps onto the x- and the y-axes and performs a fit of a Gaussian distribution to each projection, retrieving the widths of the fitted distributions, ${\sigma }_x$ and ${\sigma }_y$. Previous measurements showed that the scattering angles of individual particles along two orthogonal axes are uncorrelated, while the means of the distributions strongly correlate [17]. This suggests that ${\sigma }_x$ and ${\sigma }_y$ represent the same quantity and hence the measurements of the beam size along $x$ and $y$ can be treated as two independent measurements in order to reduce uncertainties. This procedure is mathematically equivalent to the generation of two images from the individual widths and averaging the image content. In the following, the beam size per frame is determined as the mean of the standard deviations retrieved from both fits, ${\sigma }_b=\left({\sigma }_x+{\sigma }_y\right)/2$.

4.3.2 Image generation

The beam size per frame, as the measured quantity of interest, is attributed to the positions of the linear motion stages and, in case of tomographic measurements, the rotation stage. For this, the current stage positions are queried at a 10 Hz frequency and stored with their corresponding UNIX time stamps. As this process is not synchronized to the DAQ system, a synchronisation has to be performed in the post-processing and is implemented as follows: In a first step, the position of an absorber with a sharp edge, located next to the sample, is determined. This is achieved by means of a slow, manually controlled motion of the stage, in which the beam samples the edge of the absorber, and the simultaneous observation of a change in the charge maps recorded by the detector. Subsequently, before the beginning of each two- or three-dimensional scan, the absorber is moved into and out of the beam path as a part of the automated scan sequence. With the position of the absorber edge known, the exact time of the edge transition can be located both in the data stream of stage positions and in the beam profile data, as the latter exhibits a drastic increase in deposited charge at the time of this transition. The two data streams are then correlated using their individual time stamps corrected by the determined offset. No significant drift between the time stamps has been observed. This synchronization strategy requires stable beam conditions in terms of intensity and beam position, which was measured to be satisfied with relative intensity fluctuations in the order 5% and variations of the beam position of less than 20 μm at the position of the detection plane.

With the synchronized data streams of stage positions and beam profile data, each reconstructed beam size value can be attributed to a certain point in the $x$-$y$ (2D), $x$-$y$-$\phi$ (3D) or $x$-$\phi$ (Single Slice) scan range. Equally spaced image cells within the scan range can be defined such, that either a single value is obtained per image cell, or the average of multiple values falling into a cell’s range is determined. The former enables a better image resolution, while the latter has the potential to improve the image contrast.

In case of three-dimensional measurements, the beam size obtained by a measurement as discussed above is corrected by subtracting a background value, determined as the mean observed beam width within a region where the beam did not pass through the sample. This correction compensates for the finite width of the beam after traversing the beam window and takes into account the effect of multiple Coulomb scattering in air.

Subsequently, the data set is structured in sinograms, which are representations of the corrected beam size as a function of the $x$- and $\phi$-positions for individual steps along the $y$-axis, using the assignment of beam size values to the $x$-$y$-$\phi$ phase space of the motion system. From these sinograms, slices of the material through the $x$-$z$-plane are obtained via inverse Radon transforms [30], computed individually for every step along $y$. For the measurements presented here, the open source software framework scikit-image [31]; [32] was used to perform a filtered back-projection applying a ramp filter for the reconstruction of individual slices. Subsequently, a wavelet denoising algorithm included in the framework following [33] is applied on the final images for an improvement of the image quality.

4.4 Beam characterization

The lateral size of the electron beam generated by the ARES linear accelerator was characterized in terms of a measurement as a function of the longitudinal distance from the beam exit window. For this measurement, the Timepix3 detector assembly, read out via the data acquisition system described above, was mounted to the three-dimensional $x$-$y$-$z$ linear motion stage system and a scan along the beam axis was performed covering a distance range from the beam exit window of 27 mm $<z<$ 317 mm with a step size of 10 mm. The minimum distance was limited by mechanical constraints and the maximum distance given by the range of the motion stage. Per step, the stage positions were kept constant for 60 s.

The beam size was evaluated on a per-bunch basis and the mean beam size per step is displayed in Figure 6, with the error bars representing the RMS of the beam sizes per step.

Figure 6

Figure 6. Measurement of the beam width ${\sigma }_b$ as a function of the detector distance from the beam exit window. A fit to the data represents a quadratic superposition of a focusing and two scattering contributions from the beam exit window and air.

The beam features a focal point at a distance of about 71 mm from the beam exit window with a beam size of 322 µm as a result of focusing the beam onto the sample. The trend as a function of the longitudinal distance can qualitatively be described via three terms added quadratically: the contributions represent a focusing of the beam via the quadrupoles $(f(z))$, the scattering at the beam window $(g(z))$ (both polynomials of first degree) and a contribution describing scattering in air $(h(z))$. The latter applies Equation 1 for computing the RMS lateral displacement after Equation (34.20) of [34]. This is indicated by a fit to the function

$h\left(z\right)=\sqrt{f{\left(z\right)}^2+g{\left(z\right)}^2+h{\left(z\right)}^2}=\sqrt{{\left(b_f+m_{f}x\right)}^2+{\left(b_g+m_{g}x\right)}^2+\frac{a}{\sqrt{3}}z{\theta }_0\left(z\right)}$

with the free parameters $b_{f,g}$, $m_{f,g}$ and $a$. The individual contributions are displayed in black, green and cyan. A quantitative analysis was not performed due to systematic effects such as a saturation of the detector for small beam sizes and the limited data range.

The data shows that the beam is focused at the position of the sample and the rotation axis, located at a distance of 68 mm from the beam exit window, which enables an optimal image resolution. The detection plane was positioned as close as possible, at a position of 134 mm with a beam size of 360 µm.

5 Results

5.1 Projectional measurements

In the following results, the beam size is displayed as a function of the sample position. This beam size serves as an estimator for the traversed material budget, hence a two-dimensional measurement qualitatively represents a two-dimensional projection of the material budget onto the image plane.

Figure 7 shows the result of a measurement performed in the frame-based readout mode with a scan velocity of 0.5 mm/s along the $x$-axis and 520 steps of 0.1 mm each along the $y$-axis. The measurement time amounted to 620 min.

Figure 7

Figure 7. Two-dimensional electronCT measurements of a mouse phantom, evaluated with image cell sizes of 0.05 mm $\times$ 0.1 mm (left) and 0.5 mm $\times$ 0.5 mm (right). Orange boxes indicate the regions used for the determination of the image contrast.

The data were evaluated with image cell sizes of 0.05 mm $\times$ 0.1 mm (left) and 0.5 mm $\times$ 0.5 mm (right). Both images reveal many details of the sample: the skeleton can clearly be distinguished from the tissue and is resolved to a good level of detail, exhibiting the spine, the ribs, the arms and the skull. The tissue can be well distinguished from the background. The smallest features resolved in these images are the ribs with a thickness of about 0.5 mm, which provides an upper limit for the achieved resolution.

It is readily visible, that while achieving a better resolution in Figure 7 (left), Figure 7 (right) exhibits less noise and thus a higher contrast. The contrast-to-noise ratios (CNR) are determined by

$\mathrm{C}\mathrm{N}\mathrm{R}=\frac{{\mu }_{\text{sig}}-{\mu }_{\text{bg}}}{\sqrt{{\sigma }_{\text{sig}}^2+{\sigma }_{\text{bg}}^2}}(2)$

with $\mu$ and $\sigma$ as the mean and standard deviation of values defined in a signal and a noise area. As a signal region, a homogeneous region in the shoulder of the sample representing tissue, and as a noise region, an equally sized region outside the sample have been used and are indicated as orange boxes. The CNRs result in 18.8 for the left image in Figure 7 and 34.7 for the right image in Figure 7.

5.2 Tomographic measurements

Two types of tomographic measurements have been performed using the presented phantom:

1. $x$-$y$-$\phi$ scan for a three-dimensional tomographic reconstruction.

2. $x$-$\phi$ scan for the tomographic reconstruction of a single slice in the $x$-$z$-plane.

For measurement 1, a scan velocity of 1.5 mm/s along the $x$-axis and 50 steps of 1 mm along the $y$-axis have been selected, with 37 projections recorded for a half turn, resulting in angular steps of 5°. The measurements were performed using the Katherine readout system in the data-driven readout mode, recording a continuous data stream for around 18 h.

From the three-dimensional scan, two-dimensional projection images can be generated for different rotation angles of the sample. Figure 8 shows four such projections for rotation angles of $\phi =${ 0°, 60°, 120°, 175° } with image cell sizes of 0.15 mm $\times$ 1 mm. Due to the higher scan velocity and larger step size along $y$ with respect to the images presented in Figure 7, the resolution is inferior to these. The projections acquired at rotation angles of 0° (left) and 175° (right) show similar features, but a reduced contrast can be observed in the latter. The reason for this is the positioning of the sample slightly off the rotation axis and thus the sample having different distances from the detector, which comes with an impact on the measured beam size.

Figure 8

Figure 8. Two-dimensional projections of the phantom recorded for a tomographic measurement at rotation angles of $\phi =$ { 0°, 60°, 120°, 175° } (f.l.t.r.). Orange ticks on the left side of the images indicate the vertical positions for which reconstructions are presented.

Sinograms are generated for each horizontal line, hence for every scan step along the $y$-axis, combining the data of all rotation steps. An example is shown in Figure 9, displaying the beam size as a function of $x$ and $\phi$ for a vertical position inside the head region of the phantom.

Figure 9

Figure 9. Sinogram generated from two-dimensional projections, such as shown in Figure 8 for a vertical position representing the head of the phantom.

Inverse Radon transforms are performed separately for each position along the $y$-axis, resulting in horizontal cuts, or slices, through the sample. Figure 10 shows four such slices for the vertical positions indicated via orange ticks on the left hand sides of the images shown in Figure 8. They comprise regions of the abdomen including the spine (left), the lung including the spine and a cavity for an insertable heart (center-left), the upper arm and the shoulder (center-right) and the head including the skull and the arms (right). The latter corresponds to the sinogram shown in Figure 9. Despite a larger noise contribution compared to the two-dimensional scans, the features of the phantom can be recognized and a differentiation between air, tissue-like and bone-like material is possible.

Figure 10

Figure 10. Three-dimensional electronCT measurement of a mouse phantom with an image cell size of 0.15 mm $\times$ 0.15 mm, evaluated at four positions along the vertical axis representing different anatomical regions of the phantom.

Measurement 2 was performed at a lower scan velocity of 0.7 mm/s along the $x$-axis, with 101 rotations covering a half turn. This measurement is performed at a constant $y$-position within the head (cf. Figures 9, 10 (right)) and thus yields a single horizontal cut through the sample. The measurement time was 101 min.

The result is presented in Figure 11 (right) and compared to a reconstruction slice at the same position from measurement 1 (left), both with image cell sizes of 0.15 mm $\times$ 0.15 mm. A clear improvement from (left) to (right) is visible in terms of contrast and resolution, as well as a reduction of artefacts. The CNRs, calculated from signal and background regions indicated via orange boxes and determined via Equation 2, amount to 1.96 (Figure 11, left) and 4.62 (right). These improvements are achieved by an increase in statistics: while Figure 11 (left) contains the information of 10,360 frames, (right) contains information of 57,700 frames due to the reduced stage velocity and the increased number of projections. Differences caused by a variation in the scanning procedure are neither expected nor observed.

Figure 11

Figure 11. Comparison of two tomographic reconstructions of the sample, both acquired and evaluated at the same vertical position and displayed with image cell sizes of 0.15 mm $\times$ 0.15 mm (left) shows the enlarged result of an $x$-$y$-$\phi$ scan shown in Figure 10 (right) shows the result of an $x$-$\phi$-scan with increased amount of data. Orange boxes indicate regions used for the determination of the image contrast.

It should be noted that the reconstruction via an inverse Radon transform assumes a linear dependency of the measured observable, here the beam width, on the quantity to be reconstructed, here the material budget density. While Equation 1 suggests a dependency of the opening angle and thus the beam width on the square root of the material budget, it has been found that using the detected beam width with a background subtraction as an input to the filtered backprojection, yields reasonable results with minor artefacts. Further corrections mitigating non-linearities are expected to improve the image quality and are subject to current studies.

6 Potential and limitations

The presented measurements demonstrate the technological feasibility of the electronCT imaging technique and show reasonable resolution and contrast for macroscopic objects of sizes in the order of a few tens of millimeters, resolving details such as the skeleton of a mouse. In the following, the potential and the limitations of this technology are discussed.

6.1 Scanning strategy

For the measurements presented above, a constant transverse position was chosen for the electron beam, while the phantom was moved across the transverse plane. For the studies on phantoms such a strategy is applicable whereas it would be unfeasible for imaging in medical scenarios.

Instead, dipole magnets used for steering the transverse beam trajectory could be used to move the beam relatively to a static sample. This however comes with the necessity to either move the detector synchronously with the center of the electron beam or to use a much larger detector that would be required to cover the full scan range of incident beam positions, and with the requirement of a small energy spread to avoid dispersion effects.

6.2 Spatial resolution

The achievable resolution for the electronCT technique depends on several parameters:

$•$ Step size: naturally the resolution of an image depends on the pitch of the individual image cells, which for this technique is limited by the distance between the impact positions of consecutive bunches or the step size along discretely scanned axes. This parameter also has a direct impact on the measurement time.

$•$ Beam size: the transverse size of the beam can pose a limit on the spatial resolution of the obtained image, as individual electrons are scattered at different positions at the sample. This property is limited by the accelerator performance and the geometry of the experiment.

$•$ Beam widening: already while traversing the sample, the beam size increases due to multiple Coulomb scattering. Hence, the spatial resolution suffers in the case of samples with a large thickness and high material budget.

$•$ Beam stability: shot-to-shot fluctuations or long-term drifts of the beam position and intensity at the sample can affect the spatial resolution but could partially be recovered by correcting individual image points for the center of the detected beam profile.

For the measurements shown above, the limits on the spatial resolution were posed by the RMS beam size of around 320 µm for some of the measurements, while others were dominated by the step size. The phantom was chosen sufficiently small to not expect dominating effects from a beam widening inside the sample and the position of the beam center was found to be stable within less than 20 µm.

6.3 Measurement time

A great challenge for this imaging technique lies in the measurement time. The measurements presented above range from 101 min to 18 h and are thus unacceptable for medical imaging of living beings. The measurement time however strongly depends on, and is in these studies limited by, the repetition rate of the accelerator in use. Using higher repetition rates directly decreases the measurement time, such that images such as the ones shown in Figure 7 with a repetition rate of e.g., 1 kHz could be performed within less than 10 min. For the detector in use, the readout bandwidth is limited to 85 MHits/s, such that with less than 2,000 hit pixels per frame (see Figure 4), even repetition rates of up to 42 kHz and thus measurement times of less than a minute would be theoretically possible. However, at very high rates, the synchronisation of the data acquisition with the scan of the beam position across the sample via dipole magnets requires a high precision for a proper image reconstruction.

At the same time, the requirement on the spatial resolution of an image depends on its purpose and could thus be much looser than the resolution achieved with the above studies. As the resolution impacts the number of data points and thus the total number of electron bunches required for an image, imaging could be performed much faster by adapting the target resolution.

7 Conclusion and outlook

The feasibility of the imaging technique electronCT has been demonstrated by means of two- and three-dimensional measurements. The method is based on the determination of the beam profile of a low-charge electron beam with energies in the range of few hundreds of MeV after the traversal of a phantom. Projectional images as well as tomographic reconstructions of a medical mouse phantom have been acquired using the ARES linear accelerator for beam generation and a Timepix3 detector assembly as the detection layer. The images exhibit many details of the phantom and show good resolution and contrast. This proof of concept enables studies toward an application in the context of radiation treatment with VHEE, where this method could create synergies in applying the same accelerating structure for treatment and imaging.

The presented studies expose limitations and technical challenges of the technique in the prospect of medical imaging, which lie in the measurement time and with it artifacts arising from the motion of living beings, and the spatial resolution when applied to larger objects. Hence, further studies on this technique are required to gain an understanding on the improvement potential of the measurement time, but also for an estimation of the dose delivered to a sample or patient in various imaging scenarios to explore its potential as imaging modality in a medical context. In summary, electronCT represents a candidate for in situ imaging in the context of VHEE radiotherapy and could contribute to IGRT, among others for the patient and tumor localization, under the premise of overcoming and gaining a further understanding of the aforementioned challenges.

Additional future studies serve the optimisation of this technique and comprise detailed simulations using the semiconductor detector simulation framework ${\mathrm{A}\mathrm{l}\mathrm{l}\mathrm{p}\mathrm{i}\mathrm{x}}^2$, which is capable of simulating the effects of multiple Coulomb scattering in a phantom as well as the detector response to the particle beam. Such simulations can serve as a guidance for an optimized measurement setup and strategy and provide insight into resolution limits and potential artefacts. These limits will be studied via the imaging of more generic, geometric samples such as a Derenzo phantom [35].

As indicated above, the linearity of the input to the reconstruction on the material budget is essential for an artefact-free three-dimensional imaging. To overcome this, calibration measurements applying materials of different, known properties in terms of radiation length and thickness, can furthermore be performed and applied.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

PS: Conceptualization, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Visualization, Writing–original draft. AA: Writing–review and editing. FB: Methodology, Resources, Writing–review and editing. Ld: Formal Analysis, Investigation, Writing–review and editing. HD: Investigation, Software, Writing–review and editing. KD: Formal Analysis, Investigation, Writing–review and editing. AH: Investigation, Writing–review and editing. SJ-M: Investigation, Writing–review and editing. MK: Investigation, Software, Writing–review and editing. WK: Investigation, Writing–review and editing. FM: Investigation, Software, Writing–review and editing. SR: Investigation, Writing–review and editing. SS: Conceptualization, Formal Analysis, Investigation, Methodology, Software, Writing–review and editing. TV: Investigation, Writing–review and editing. HW: Methodology, Supervision, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This project received funding via the DESY Generator Program.

Acknowledgments

The phantom has been obtained from Universitätsklinikum Hamburg-Eppendorf. The authors acknowledge support from DESY (Hamburg, Germany), a member of the Helmholtz Association HGF, and thank all the technical groups at DESY for their work and support in the ARES implementation, maintenance and operation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Ronga MG, Cavallone M, Patriarca A, Leite AM, Loap P, Favaudon V, et al. Back to the future: very high-energy electrons (vhees) and their potential application in radiation therapy. Cancers (2021) 13:4942. doi:10.3390/cancers13194942

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Farr JB, Parodi K, Carlson DJ. Flash: current status and the transition to clinical use. Med Phys (2022) 49:1972–3. doi:10.1002/mp.15401

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Vozenin MC, Bourhis J, Durante M. Towards clinical translation of flash radiotherapy. Nat Rev Clin Oncol (2022) 19:791–803. doi:10.1038/s41571-022-00697-z

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Lin B, Gao F, Yang Y, Wu D, Zhang Y, Feng G, et al. Flash radiotherapy: history and future. Front Oncol (2021) 11. doi:10.3389/fonc.2021.644400

PubMed Abstract | CrossRef Full Text | Google Scholar

5. De Los Santos J, Popple R, Agazaryan N, Bayouth JE, Bissonnette JP, Bucci MK, et al. Image guided radiation therapy (igrt) technologies for radiation therapy localization and delivery. Int J Radiat Oncology*Biology*Physics (2013) 87:33–45. doi:10.1016/j.ijrobp.2013.02.021

CrossRef Full Text | Google Scholar

6. Ruchala KJ, Olivera GH, Schloesser EA, Mackie TR. Megavoltage CT on a tomotherapy system. Phys Med and Biol (1999) 44:2597–621. doi:10.1088/0031-9155/44/10/316

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Molière G. Theorie der Streuung schneller geladener Teilchen II - Mehrfach-und Vielfachstreuung. Z Naturforsch Teil (1948) A(3):78–97. doi:10.1515/zna-1948-0203

CrossRef Full Text | Google Scholar

8. Bethe HA. Molière’s theory of multiple scattering. Phys Rev (1953) 89:1256–66. doi:10.1103/PhysRev.89.1256

CrossRef Full Text | Google Scholar

9. Highland VL. Some practical remarks on multiple scattering. Nucl Instrum Meth (1975) 129:497–9. doi:10.1016/0029-554X(75)90743-0

CrossRef Full Text | Google Scholar

10. Lynch GR, Dahl OI. Approximations to multiple coulomb scattering. Nucl Instrum Meth B (1991) 58:6–10. doi:10.1016/0168-583X(91)95671-Y

CrossRef Full Text | Google Scholar

11. Lutz G. Semiconductor radiation detectors. Verlag Berlin Heidelberg: Springer (1999).

Google Scholar

12. Spannagel S, Wolters K, Hynds D, Tehrani NA, Benoit M, Dannheim D, et al. Allpix²: A modular simulation framework for silicon detectors. Nucl Instrum Meth A (2018) 901:164–72. doi:10.1016/j.nima.2018.06.020

CrossRef Full Text | Google Scholar

13. Agostinelli S, Allison J, Amako K, Apostolakis J, Araujo H, Arce P, et al. GEANT4 - a simulation toolkit. Nucl Instrum Meth A (2003) 506:250–303. doi:10.1016/S0168-9002(03)01368-8

CrossRef Full Text | Google Scholar

14. Spannagel S, Wennlöf H, Schütze P. Allpix squared - generic pixel detector simulation framework (2023). doi:10.5281/zenodo.8171803

CrossRef Full Text | Google Scholar

15. Weinrich U, Fuchs R, Sust E. Assembly of the carbon beam gantry at the heidelberg ion therapy (HIT) accelerator. Conf Proc C (2008).

Google Scholar

16. Chinniah S, Deisher AJ, Herman MG, Johnson JE, Mahajan A, Foote RL. Rotating gantries provide individualized beam arrangements for charged particle therapy. Cancers (2023) 15:2044. doi:10.3390/cancers15072044

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Jansen H, Schütze P. Feasibility of track-based multiple scattering tomography. Appl Phys Lett (2018) 112:144101. doi:10.1063/1.5005503

CrossRef Full Text | Google Scholar

18. Burkart F, Aßmann R, Dinter H, Jaster-Merz S, Kuropka W, Mayet F, et al. The ARES linac at DESY. In: Proc. LINAC’22 International Linear Accelerator Conference; Liverpool, UK. Geneva, Switzerland: JACoW Publishing (2022). p. 691–4.

Google Scholar

19. Poikela T, Plosila J, Westerlund T, Campbell M, Gaspari MD, Llopart X, et al. Timepix3: a 65k channel hybrid pixel readout chip with simultaneous toa/tot and sparse readout. J Instrumentation (2014) 9:C05013. doi:10.1088/1748-0221/9/05/C05013

CrossRef Full Text | Google Scholar

20. Bienvenu G. Dark current in high-gradient accelerator. IEEE Trans Electr Insul (1989) 24:1037–9. doi:10.1109/14.46334

CrossRef Full Text | Google Scholar

21. Dudak J. High-resolution x-ray imaging applications of hybrid-pixel photon counting detectors timepix. Radiat Measurements (2020) 137:106409. doi:10.1016/j.radmeas.2020.106409

CrossRef Full Text | Google Scholar

22. Turecek D, Jakubek J, Trojanova E, Sefc L, Kolarova V. Application of timepix3 based cdte spectral sensitive photon counting detector for pet imaging. Nucl Instr Methods Phys Res Section A: Acc Spectrometers, Detectors Associated Equipment (2018) 895:84–9. doi:10.1016/j.nima.2018.04.007

CrossRef Full Text | Google Scholar

23. Pitters FM, Alipour Tehrani N, Dannheim D, Fiergolski A, Hynds D, Klempt W, et al. Time and energy calibration of Timepix3 assemblies with thin silicon sensors. CLICdp-Note-2018-008 (2018).

Google Scholar

24. Burian P, Broulím P, Jára M, Georgiev V, Bergmann B. Katherine: ethernet embedded readout interface for timepix3. J Instrumentation (2017) 12:C11001. doi:10.1088/1748-0221/12/11/C11001

CrossRef Full Text | Google Scholar

25. Mánek P, Burian P, David-Bosne E, Smolyanskiy P, Bergmann B. Track Lab: extensible data acquisition software for fast pixel detectors, online analysis and automation (2023).

Google Scholar

26. Wegner M, Frenzel T, Krause D, Gargioni E. Development and characterization of modular mouse phantoms for end-to-end testing and training in radiobiology experiments. Phys Med and Biol (2023) 68:085009. doi:10.1088/1361-6560/acc566

CrossRef Full Text | Google Scholar

27. Thorlabs Inc. 300 mm linear translation stage with integrated controller, stepper motor (2024). Available from: https://www.thorlabs.com/newgrouppage9.cfm?objectgroup_id=7652&pn=LTS300/M (Accessed June 10, 2024).

Google Scholar

28. Thorlabs Inc. PRM1/MZ8 - Ø1 motorized precision rotation stage (metric) (2024). Available from: https://www.thorlabs.com/thorproduct.cfm?partnumber=PRM1/MZ8 (Accessed June 10, 2024).

Google Scholar

29. Dannheim D, Dort K, Huth L, Hynds D, Kremastiotis I, Kröger J, et al. Corryvreckan: a modular 4d track reconstruction and analysis software for test beam data. J Instrumentation (2021) 16:P03008. doi:10.1088/1748-0221/16/03/P03008

CrossRef Full Text | Google Scholar

30. Deans S. The Radon transform and some of its applications. In: Dover books on mathematics series. Dover Publications (2007).

Google Scholar

31. scikit-image development team. scikit-image (2024). Available from: http://scikit-image.org/ (Accessed June 10, 2024).

Google Scholar

32. van der Walt S, Schönberger JL, Nunez-Iglesias J, Boulogne F, Warner JD, Yager N, et al. scikit-image: image processing in Python. PeerJ (2014) 2:e453. doi:10.7717/peerj.453

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Chang S, Yu B, Vetterli M. Adaptive wavelet thresholding for image denoising and compression. IEEE Trans Image Process (2000) 9:1532–46. doi:10.1109/83.862633

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Workman RL, Burkert VD, Crede V, Klempt E, Thoma U, Tiator L, et al. Review of particle physics. PTEP 2022 (2022) 2022:083C01. doi:10.1093/ptep/ptac097

CrossRef Full Text | Google Scholar

35. Derenzo SE, Budinger TF, Cahoon JL, Huesman RH, Jackson HG. High resolution computed tomography of positron emitters. IEEE Trans Nucl Sci (1977) 24:544–58. doi:10.1109/tns.1977.4328738

CrossRef Full Text | Google Scholar