This invention relates to a method and an apparatus for proton computed tomography, and in particular, methods for the 3-dimensional reconstruction in computed tomography (CT). It has specific, but not exclusive, application to computerized tomographic image capture and reconstruction for applications in proton and other charged particle therapies in the treatment of cancer.

The most frequently used charged particle therapy uses protons, so the description that follows will refer exclusively to protons. However, other forms of charged particle or hadron radiation imaging and therapy can be used. Embodiments of this invention may make use of other types of radiation.

A potential patient for proton radiotherapy would typically receives an x-ray CT scan from which the subsequent treatment would be planned. Planning requires the accurate estimation of proton range (a function of its energy) within the patient. However, x-rays and protons interact with tissue through different physical mechanisms. The current uncertainty in proton range as compared with x-rays is commonly estimated as ±3.5% (at 95% confidence limit), so if the beam passes through 20 cm of tissue, then the Bragg Peak could be anywhere within ±7 mm. Using protons instead of photons (x-rays) to produce CT imagery greatly reduce these range errors.

There is, in addition, a need for on-treatment imaging to assist in ensuring that patients are appropriately positioned during treatment and to accommodate changes in the internal anatomy of the patient during a course of treatment. Although CT using proton beams (here referred to as proton-CT, pCT) has a long history, it has recently received revived interest due to the need to increase conformance of radiotherapeutic treatment by reducing delivery errors.

The major effort has been on the reconstruction of an image from the relative stopping power in a patient as this has the potential to reduce proton range uncertainties and therefore be more accurate and offer reduced margins in treatment planning. Early work on pCT (Cormack A M, Koehler A M, 1976, Quantitative proton tomography: preliminary experiments, Phys Med Biol, 21, 560-9), took a similar approach as conventional x-ray CT, integrating the proton flux after it had passed through a patient (or in experimental work, a surrogate object - a "phantom"). However, protons experience multiple Coulomb scattering (MCS) as they pass through material so the resulting spatial resolution in the reconstructed pCT is poor and not of sufficient quality for clinical application. Later work has concentrated on proton-tracking approaches where the trajectories of individual protons are determined before and after the patient.

The general arrangement of a proton-tracking pCT system is shown in Figure 1, where the patient 1 is situated between two pairs of position-sensitive detectors (PSD) 3, 4; 5, 6; together with a residual energy-range detector (RERD) 7. There is a relative rotation between the patient and apparatus about the system's isocentre 2. The trajectory of a single proton is shown at 8, with an illustrative non-linear path through the patient. The input and output trajectories are recorded (where 'X' marks the measured position in each PSD 3, 4; 5, 6). The energy of the incident protons needs to be increased beyond the normal treatment energy to ensure that a majority of incident protons will pass through the patient. The RERD 7 records the residual energy of the individual tracked protons. The PSDs 3, 4; 5, 6 can include any segmented sensor, such as crossed plane of silicon strip sensors, multi-wire proportional chambers or scintillating-fibre hodoscopes, capable of recording, with sufficient precision, the x-y coordinates of a transient proton. The RERD 7 can be based on several technologies such as crystal scintillator calorimeters, plastic scintillators stacks or layers of silicon sensors (strips or active pixel sensors). To calculate the proton energy lost within the patient, ideally there would be a detector that records the energy of each incident proton. This is not currently possible so incident energies are estimated from the properties of the proton beam and the generating source. To reconstruct a pCT image, the x-y coordinates recorded by the four PSDs 3, 4; 5, 6, namely for the i ^th proton, together with its recorded residual energy, (ER)I are used. These quantities are indicated in Figure 2. The entry and exit locations at the patient's surface, respectively are estimated by extrapolating the two straight lines derived from the PSD coordinates and a knowledge of the geometry of the system. The energy (ΕΑ)i absorbed in the patient for the i ^th proton is derived from estimating its incident energy minus its residual energy, (ER)i The process is repeated many thousands of times for a fixed orientation of the patient with respect to the apparatus. For a full volumetric pCT reconstruction, the relative orientation about the isocentre is changed by one or a few degrees and the process of recording repeated. The process, in turn, is repeated until the total relative movement equates to typically 180°, but could be less or up to 360°.

The interaction of energetic protons with matter is chiefly through Coulomb interactions which is a nuclear process. Secondly, inner-shell as well as outer-shell electrons contribute to energy-loss except at very low proton kinetic energies. The protons lose most of their energy near the end of their paths - the so-called Bragg Peak. Other advantages of charged particles, such as protons, are their reduced lateral dose penumbra as well as absence of an exit dose tail in the tissue, which is after the Bragg Peak. Proton energy losses, which are statistical processes, gradually bring the particle to rest. The fluctuation in the proton range is termed "range straggling". The other significant parameter in radiotherapy is the lateral broadening of incident beam as particle penetrates into the patient. The multiple small-angle deviations along the proton path, in total, lead to a change in particle directions that result in a deviation in primary path. The lateral broadening of beam due to multiple scattering is often called lateral scattering. In addition, a fraction of proton experience non-elastic nuclear interactions with the target nuclei, so the fluence of the beam reduces with depth into the patient. The basis of this invention is that these phenomena can be exploited to yield alternate forms of pCT.

The primary aim on this invention is to provide a benefit in the practical application of proton therapy. The major focus has been on the reconstruction of relative stopping power in a patient to reduce proton range uncertainties and therefore margins in treatment planning. There are also potential advantages for guiding and modulating proton beams during treatment: often termed image-guided therapy.

In addition, embodiments of the invention pCT can yield complementary information. For example, volumetric images of relative stopping power, scattering-power and attenuating power in the patient may all be desirable for accurate treatment planning.

From a first aspect, this invention provides method of 3-dimensional image reconstruction in computerised tomography comprising: a. directing a plurality of particles at a target object from a plurality of incident angles; b. for each incident particle, recording its position and trajectory before and after the target object; c. employing a computerised tomography reconstruction algorithm to transform the recorded data for each incident angle into a backprojected computerised tomography image; wherein d. the computerised tomography image reconstruction algorithm is based on one of or a combination of two or more of the stopping-power of the particles, the scattering-power of the particles, the attenuating-power of the particles or straggling-power of the particles; and e. the computerized tomography image reconstruction algorithm converts the power or powers of step d. to a water-equivalent path length (WEPL) for each incident angle after the backprojection of that projection.

The advantages in exploiting these different quantities to produce pCT images is that the resulting images possess identifiable improvements in some image quality metrics (for example, improved image contrast) and, more importantly, can require much reduced instrument complexity and cost for their measurement compared to the conventional stopping power pCT.

The presence in a patient's body of many different tissues, each with its own density and chemical composition, constitutes a challenge in radiotherapy, since reconstruction of an image would require simulating the radiation's progression through each of these materials. The usual approach to this problem is to apply the concept of water equivalent path length (WEPL), defined as the radiological depth between a source and any other point of interest or calculation point, as determined by the linear attenuation of each material in the path. In other words, if we consider a single particle traversing several tissues of different thicknesses and densities, the WEPL concept scales all these tissues to the depth of water that has the same attenuating effect. For protons, WEPL is the thickness of water that would produce an equivalent energy-loss, scattering, attenuation or straggling. A tomographic acquisition determines a set of WEPL data (projection data). The relationship between a given kinematic property of the protons and the estimated WEPL must be known. This relationship can be determined for differing materials by such methods as Monte Carlo simulations or experimental measurements.

Normally, in producing a CT image, the measured property is converted to its WEPL equivalent prior to applying the backprojection algorithm. Here, it is proposed here that this conversion is applied after backprojection.

Methods embodying the invention for reconstructing scattering-power and attenuating-power have no requirement for any knowledge of the exit energy of the protons. This means that pCT for these two quantities is possible with only trackers (PSDs) and no range-telescope or calorimeter (i.e. no RERD). The absence of a RERD enables less complex, less expensive and more compact pCT systems.

In step d., the algorithm may be further based on the incident energy of the particles. The incident energy of the particles is detected or it may be estimated.

The reconstructed image is a 2-dimensional slice image or it may be a full 3- dimensional volumetric image.

The particles are typically protons, but may be other hadrons or other charged particles.

The trajectories of individual particles before and/or after traversal of the target object may be recorded by a multiplicity of position-sensitive detectors that are encountered by the particles before they encounter the object.

In step b, one or more of the kinetic energy of individual particles after traversal of the target object, the multiple Coulomb scattering (MCS) angle of individual particles after traversal of the target object may be recorded, the fraction of individual particles that transverse the target object and the deviations of the energy of a particle after it has transversed the target object may additionally be recorded.

The method typically requires backprojecting the or each kinematic variable of interest (power) rather than WEPL, and only converting to WEPL after backprojection so that conversion to WEPL does not happen after the generation of the final backprojection image but is applied to each angle after back projection of that projection.

A method of image reconstruction according to any preceding claim in which the backprojected image may be expressed as: where M is a set of projections recorded in step b, N _m is a set of particles in the mth projection and Δβ _Μ is the angular increment associated with the projection

The backprojected image in terms of stopping-power may be expressed as:

where is the energy loss for the nth proton.

The backprojected image in terms of scattering-power may be expressed as:

where is the scattering angle for the nth proton.

The backprojected image in terms of attenuating-power may be expressed as: where L _m is the set of all protons absorbed in the object and N _m is the set of all protons that were not absorbed in the target object.

The backprojected image in terms of straggling-power may be expressed as:

Production of a reconstructed image in terms of WEPLs is achieved through the function, where the WEPL for a specified kinematic property, % is given by:

where is the change in the kinematic property of protons

is backprojected for each angle and then converted to WEPL, prior to summing the

various projection contributions. Though not essential for reconstructing stopping- power, it is important for the other powers.

A method of image reconstruction according to any preceding claim in which the target object may be a living organism such as a living person.

A method of 3-dimensional image reconstruction in computerised tomography substantially as hereinbefore described with reference to Figures 3 to 6 of the accompanying drawings.

Computerised tomography apparatus comprising a source of particles, first detection means for detecting the position and energy of particles prior to their traversing an object, second detection means for detecting the position and energy of particles prior to their traversing an object, and analysing means configured to reconstruct an image of the object by analysis of output of the detection means using a method according to any preceding claim.

Embodiments of the invention will now be described in detail, by way of example, and with reference to the accompanying drawings, in which:

Figure 1 shows diagrammatically the general arrangement of a proton-tracking pCT system;

Figure 2 shows, in more detail, such instrumentation with the notation used for measured quantities and the employed reference axes notation; Figure 3 shows diagrammatically the trackers-only instrumentation required for pCT reconstruction based on scattering-power and attenuating-power being an embodiment of the invention;

Figure 4 shows the simulated pCT reconstructions of a phantom based on the four different pCT, In the order relative, stopping, scattering, attenuating and straggling power, described below;

Figure 5 shows simulated pCT reconstructions for stopping and straggling power with low numbers of protons; and

Figure 6 shows experimental images of relative scattering-power in two axial slices through a 75 mm diameter PMMA phantom.

A system for implementing methods embodying the invention is shown in Figure 3. The system similar to that of Figure 1, described above, but with the residual energy- range detector (RERD) 7 omitted.

The theory underlying embodiments of the invention will first be described. The behaviour of an individual proton in its interaction with matter is not deterministic, in that its state is not fully governed by some set of parameters and initial conditions. However, an ensemble of protons will behave in a reproducible, quasi- deterministic manner. After passing through a thickness of material, protons of the same initial energy will typically have energies close to a mean value. Around this mean there will be a spread of actual values, due to energy straggling, that will increase with depth. Not only do protons lose energy as they pass through a patient, they also scatter laterally, and a proportion is lost from the beam through inelastic interactions between the protons and atoms within the patient, so there is an increasing spread of energies and loss of protons within the beam as it progresses. So there are four quantities that describe ensembles of protons in a beam that has passed through a target: stopping-power (energy-loss), scattering-power (lateral deflection), attenuating-power (loss of protons due to nuclear interactions) and straggling-power (increase in energy spread).

The desirability of using not only relative stopping-power, but also relative scattering- power and attenuating-power is known. The methods presented here enable all four quantities to be determined experimentally with high spatial-resolution. However, it is stopping-power that is the most crucial quantity to determine for proton therapy planning.

Each of relative stopping, scattering and attenuating-powers have differing dependences on the elemental composition of a material (i.e., effective atomic number,

The new methods of obtaining pCT have some additional advantages. For example, scattering-power reconstruction does exhibit greater contrast but at the price of increased noise than stopping-power CT. It is possible to combine two or more pCT images to yield an improved quality pCT. It is possible to formulate image reconstruction algorithms for all four quantities with high spatial resolution using estimated non-linear proton paths. The principle is demonstrated using Monte Carlo simulations. Note that the data set required for reconstructing scattering-power, attenuating-power and straggling-power never exceeds that needed for reconstructing stopping-power pCT. Therefore, these complementary quantities can be found without any additional radiation exposure to the patient.

In particle-resolving pCT, protons are used as a probe of material property, and the kinematic properties of individual protons are determined or estimated. For present purposes, the state of a proton, Sj, will be considered fully specified by the generalized set of kinematical properties, where the z-axis is the beam direction, x andy are the lateral dimensions, a and β are the angles with respect to the respective lateral directions and E is the energy of the proton. The state of each proton must be determined prior to entry into a target fj = in) and, if the proton is transmitted rather than absorbed or backscattered, also upon exit fj = out).

Reconstruction of an image from stopping-power and straggling-power require knowledge of the full states, Sm and Sout, but to reconstruct an image from the scattering-power and attenuating-power, however, Eout is not needed. This means that pCT of scattering and attenuation can be performed with only the pairs of position sensitive detectors (PSDs) 3,4; 5,6 as depicted in Fig. 3, as there is no requirement for a RERD.

The proton power of type, τ, at a spatial position, r, can be expressed as, where v is the direction vector of the proton and is the expectation-value of

the appropriate kinematic property. For stopping power (s), scattering (t), attenuating (i/) and straggling (v) powers, the kinematic parameters are, respectively,

Here, E is the proton kinetic energy, Θ is the scattering angle, φ o and φ are initial proton fluence and fluence at some depth, and, (E) is the expectation-value of proton kinetic energy. Here Θ is defined as the total scattering angle, but it may alternatively be defined as the projected angle.

Although a dependence on proton energy has not been explicitly stated in the arguments of the proton power in eqn. (2), some dependence will generally be present. The possibility of pCT relies on the fact that the relative proton power is only weakly dependent on proton energy.

The generalised relative power for a material can be defined as:

where refers to the respective power for water. This definition

as a relative quantity is directly analogous to the definition of Hounsfield Unit in x-ray CT, which is defined in terms of the linear attenuation coefficient relative to water.

Unlike stopping-power and attenuating-power, which can be considered as well- defined quantities dependent only on the elemental composition at a given point, scattering-power is a non-local quantity. However, relative scattering-power can be considered local.

The mechanism underlying the attenuating power is the inelastic nuclear cross-section of the constituent elements. This cross-section exhibits threshold and resonance behaviour, decaying to a plateau at sufficiently high proton energies. Relative attenuating-power can only be considered an energy-independent quantity in the plateau region, the onset of which increases with atomic number.

Typically, CT reconstruction problems are formulated in terms of the water-equivalent path-length (WEPL) through an object. For protons this is thickness of water that would produce an equivalent energy-loss, scattering, attenuation or straggling effect. Quite generally, we can define WEPL for a power, τ, as

where the integral over r' picks out the location of the proton at a given path-length and the integral over λ is a line-integral over this path.

A tomographic acquisition produces a set of WEPL data (the so-called "projection data"). Tomographic reconstruction can be thought of as a mechanism to "undo" the line-integral to provide the quantity ρ(τ), which is the CT image.

There needs to be some known relationship between a change in a kinematic property of protons, and the estimated WEPL. In general, this mapping can be expressed

as:

What is required in practice, is the inverse mapping:

where denotes the inverse function of f This relationship can be arrived at

using theoretical formulae, Monte Carlo simulations or experimental calibrations.

The process of CT reconstruction is described in terms of the Backprojection-then- filtering (BPF) algorithm proposed by Poludniowski et al (Poludniowski G, Allinson N M and Evans P M, 2014, Proton computed tomography reconstruction using a backprojection-then-filtering approach, Phys Med Biol. 59, 7905-18), though other methods are possible. The structure of the algorithm for a set of pixels *^ in a 2- dimensional axial slice, is:

where bj is the backprojected image of WEPL and kj is a kernel. The kernel expressed in the form,

where:

where J _n and H _n are nth order Bessel and Struve functions, respectively.

The projected WEPL for the nth proton, is related to the relative-power distribution in voxels, by,

where a _n j is a system matrix. The matrix elements take values equal to the path-length of the nth proton through the ; ^' th voxel. Since the ray will miss most voxels, the matrix is sparse in The backprojected image can be expressed as,

where Μ is the set of projections, N _m is the set of protons in the mth projection and is the angular increment associated with the projection

This generalised derivation can be satisfied for the above-mentioned kinematical properties of the proton beam, as discussed below.

Backprojection for stopping-power. The application of the BPF algorithm to reconstructing proton stopping-power is relatively straightforward, because the energy-loss of a proton passing through an object is quasi-deterministic. Each proton typically provides a good estimate of the expectation value of energy loss and therefore a plausible estimate of the path length travelled. The backprojection operation of equation (14) may be approximated as,

where is the energy-loss for the nth proton. Strictly, however,

The backprojection operation is itself a voxel-specific

estimator for the stochastic variable being backprojected, so this allows the formulation of the backprojection as:

That is, backproject ΔΕ _η at each angle rather than w _n and then convert to WEPL prior to summing the various projection contributions. The reordering of the conversion to WEPL that constitutes the difference between eq. (15) and (16) is important for all powers.

Backprojection for scattering power. In this case, scattering-angle is a random variable. The increase in the square of scattering angle due to passage through an object cannot be treated as a quasi-deterministic variable. It is therefore not advisable to take the

square-angle of a single proton and assume because this would typically

result in a biased estimate. A more appropriate form for the backprojection is:

Backprojection for attenuating-power. Attenuating-power can be treated in a similar manner to the cases above. However, there is the added complication of the normalising fluence, φο. In some cases, it may be sufficient to estimate values based on known incident beam properties. This assumption is not required if φο is calculated by backprojection. The backprojection relation becomes:

where L _m is the set of all protons absorbed in the patient and N _m is the set of all proton that were not absorbed in the patients. Note that the most-likely path for the proton in set Lm will be a straight-line, as insufficient data is available to estimate a curved path.

Backprojection for straggling-power. The backprojection of straggling-power has a different complication. Since, a linear combination of backprojections is necessary. The relation becomes:

To demonstrate the range of pCT methods, simulations using a well-known particle physics Monte Carlo simulation package (FLUKA: http://www.fluka.org/fluka.php) were undertaken. The test object (phantom) was a 180 mm diameter water cylinder of length 40 mm. Embedded within the phantom were four cylinders of 20 mm diameter and length, each comprising of skeletal muscle, cortical bone, adipose or air. Additionally, a tungsten bar (diameter: 1 mm; length 20 mm) was inserted into the phantom to allow quantification of spatial resolution. The source was modelled as a parallel-beam of mono-energetic protons with a kinetic energy of 200 MeV. A total of 180 projections were simulated equi-spaced over π radians with a fluence of 300 protons per mm ² in each projection.

A modified version of the "Backprojection then Filtering" (BPF) CT reconstruction was used (Poludniowski G, Allinson N M and Evans P M, 2014, Proton computed tomography reconstruction using a backprojection-then-filtering approach, Phys Med Biol. 59, 7905- 18). Instead of WEPL being backprojected, the appropriate quantity from equations (16) to (19) was used. Conversion to WEPL was carried out after the backprojection of all histories in a projection. Fig. 4 (a - d) shows the central slice for each case. Of the four material inserts, the top sample (muscle) is only discernible in the stopping-power image, due to the windowing, intrinsic contrast and relative noise levels in the images. While images (a) to (c) display similar characteristics, the enhancement at boundaries between inhomogeneities is very prominent in (d) (straggling-power).

Reconstructions of straggling-power show high sensitivity to small high contrast objects. Figure 5a shows a reconstruction of stopping- power for the phantom with the proton fluence reduced to 5 protons per mm ² per projection (corresponding to scan dose of 80 μGy at the phantom centre). The tungsten bar is no longer resolvable due to increased noise. Figure 5b shows the corresponding reconstruction for straggling- power, in which the tungsten bar is clearly visible.

In addition to simulations, Figure 6 shows experimental images of relative scattering- power in two axial slices through a 75 mm diameter PMMA phantom. The proton beam energy was 125 MeV with a diameter of 85 mm. Figure 6a shows a 1 mm diameter tungsten carbide sphere, while Figure 6b shows three material inserts.