WO2013119887A1

WO2013119887A1 - Method and statistical validation technique for detecting differences between radiation therapy images with application to the detection and control of radiation therapy treatment delivery errors

Info

Publication number: WO2013119887A1
Application number: PCT/US2013/025245
Authority: WO
Inventors: John James Gordon; Jeffrey V. Siebers
Original assignee: Virginia Commonwealth University
Priority date: 2012-02-09
Filing date: 2013-02-08
Publication date: 2013-08-15

Abstract

A method and statistical validation technique detects differences between radiation therapy verification images with application to the detection and control of radiation therapy treatment delivery errors. The process takes reference and measured images of radiation therapy treatment fields and identifies deviations between them, for purposes of detecting and correcting delivery errors that could degrade treatment effectiveness or cause injury. For each pixel or voxel, the percent dose difference between measured and reference images is computed, and scaled by a quantity that depends on the local dose gradient, the measured image noise, and the measured positional variability of field-defining elements. The measured image is classified as error-free if the percent of scaled pixel values below the threshold exceeds a cutoff. The invention further comprises a method of measuring the noise and positional variability of field-defining elements, either from a single image or a population of images, as input to the detection algorithm. The invention further comprises a method of quantifying detection performance, allowing detection thresholds and cutoffs to be optimized.

Description

METHOD AND STATISTICAL VALIDATION

TECHNIQUE FOR DETECTING DIFFERENCES BETWEEN RADIATION THERAPY IMAGES WITH APPLICATION TO THE DETECTION AND CONTROL OF RADIATION THERAPY TREATMENT DELIVERY

ERRORS

DESCRIPTION BACKGROUND OF THE INVENTION

Field of the Invention

The invention pertains generally to radiation therapy employed for curative or palliative cancer treatment. More specifically, the invention pertains to the imaging of treatment beams or fields that constitute a radiation therapy treatment plan, for purposes of verifying the correctness of such fields prior to or during treatment delivery. The invention comprises a process that takes a pair of images (e.g., reference and measured images of treatment fields) and identifies deviations between them, for purposes of detecting and correcting delivery errors that could degrade treatment effectiveness or cause injury. The invention further comprises a statistical validation technique to quantify detection accuracy, enabling calibration of the detection algorithm. Background Description

External beam radiation therapy (EB T) is one of the main treatment modalities for cancer. It generally involves delivering multiple shaped and modulated radiation beams or fields to a tumor or target location in the patient's body. The patient is positioned on a table. Radiation beams, consisting of x-rays or high energy electrons or other particles, are generated by an accelerator and directed at the target from different angles. The beams converge at the target, resulting in significant cumulative radiation in and around the target tissue, and lesser amounts of radiation to peripheral tissues traversed by the beams. Variable beam angles are achieved by directing radiation from a movable gantry, or by moving the table, or both. A common arrangement is shown in Fig. 1. In this case the gantry rotates in a vertical plane around a fixed point called the isocenter, and the table rotates in a horizontal plane around the isocenter.

More particularly, Fig 1 shows a linear accelerator used to deliver radiation therapy. Treatment fields or beams 1 1 are delivered from a gantry 12 that rotates around a horizontal axis B. The patient is supported on a table 13 that rotates around a vertical axis A. Axes A and B intersect at the fixed isocenter 1. The line from the radiation source 14 located in the gantry through the isocenter is the beam's central axis. An electronic portal imaging device (EPID) 15 is shown retracted in the base of the gantry 12. When deployed, the EPID panel extends on the far side of the patient, perpendicular to the central axis and centered on the central axis, so that it can capture dosimetric images of the delivered beam. EPID images can be captured with or without the patient present.

The patient's customized treatment plan, created using a treatment planning system (TPS), could call for radiation to be delivered via a number of fixed beams, each corresponding to a different combination of gantry and table rotations. In this case, the beam remains off while gantry and table are rotated to the next position, the beam is delivered with gantry and table stationery, then the gantry and table advance to the next position, and so on. Alternatively, the treatment plan could call for radiation to be delivered via some number of arcs, in this case, the beam remains off while the table is rotated to the next position, the beam is delivered continuously while the table remains stationery and the gantry rotates in an arc around the target, then the table advances to the next position, and so on. Alternative means of delivery exist, such as accelerators mounted on robotic arms that can translate and rotate freely. However, the principle of delivering treatments via multiple fixed beams or arcs remains the same.

Radiation dose refers to energy deposited in patient tissue by radiation. At sufficient levels, deposited dose will kill living cells, including cancerous cells and non-cancerous normal tissue cells. Different cancers and normal tissues (i.e., organs) exhibit different sensitivities to radiation. The primary goal of radiation therapy is to deliver sufficient dose to the target to kill cancerous cells (or retard their growth in the case of palliative treatment). A secondary goal is to minimize the dose delivered to surrounding normal tissues, so as to minimize the adverse consequences to the patient's healthy organs. In general, dose to normal tissues is unavoidable. The goal of treatment planning is to find a beam configuration that optimizes the tradeoff between target and normal tissue doses for the patient's specific disease distribution and anatomy. Typically this is achieved via a mathematical optimization, which finds a combination of beams that generate a dose distribution that is as conformal as possible to the target, i.e., falls off as sharply as possible outside the target. For different cancers and treatment sites, criteria exist for minimum dose to the target, and maximum tolerable doses to normal tissues.

In modern radiation therapy, conformal dose distributions are achieved by means of beam shaping and modulation. Planning techniques which embody this approach include intensity modulated radiation therapy (IMRT) (Bortfeld, T., "IMRT: a review and preview", Physics in Medicine and Biology, 51, R363-R376, 2006) and volumetric modulated arc therapy (VMAT). (Otto, ., "Volumetric modulated arc therapy; IMRT in a single gantry arc", Med. Phys., 35, 310, 2008) Fig, 2 shows a common method of beam shaping and modulation. The central axis of the accelerator is a line passing from the beam source through the isocenter. In a plane perpendicular to the central axis, the primary beam generated by the accelerator is approximately radially symmetric. This radially symmetric beam profile is measured, and entered into the treatment planning system, and becomes the basis for treatment planning.

More particularly, Fig. 2 is a cross-section through a linear accelerator treatment head, patient and EPID. The view is perpendicular to the direction of travel of the Y jaws 21 , and parallel to the direction of travel of X jaws 22 and multileaf collimator (MLC) leaves 23. The shaped beam 24 is defined by jaws and MLC leaves 23. The delivered beam, either with or without the table, patient or phantom present, can be imaged by the extended EPID. The EPID detector array 25 consists of an array of solid state detectors, which record delivered dose.

Under control of the delivery system, the primary beam is truncated to a rectangle by two sets of movable opposed jaws (X jaws and Y jaws). Further shaping of the beam is achieved by a MLC. The MLC consists of two banks of opposed leaves, which can independently move in and out of the field to create irregular beam shapes, or apertures. Additionally, both the jaws and MLC are fixed to a rotating collimator that allows them to be rotated around the central axis, in a plane perpendicular to the central axis. The jaws, MLC and collimator are all housed within the gantry. At an instant of time during delivery of a beam, the beams-eye view of the target (i.e., the view from the source to the target along the central axis) will appear as in Fig. 3. The jaws 31, 32 and MLC leaves 33 combine to create a shaped beam aperture. More particularly, Fig 3 is a beams-eye view from the linear accelerator source along the central axis, showing the circular outline of the rotatable collimator, the X and Y jaws 31 ,32 which define the outer limits of the shaped beam, and the upper and lower banks of opposed MLC leaves 33, which move to define the irregular field shape at each point during beam delivery.

Finally, under control of the delivery system, MLC leaves 33 and jaws 31,32 can move at variable speeds during beam delivery to superimpose the radiation delivered by a number of consecutive beam apertures. The end result is a modulated two-dimensional radiation profile as in Fig.4, which exhibits "peaks" and "valleys" (i.e., regions of higher and lower intensity) in the plane perpendicular to the central axis. Fig. 4 is an example of the dose produced by a modulated treatment beam. The treatment planning system will generate a number of such modulated beams, typically delivered from different gantry and table angles. Use of modulated beams permits the TPS to achieve a cumulative patient dose distribution that is more conformal to the target than would be possible with simpler, unmodulated beams.

The overall success of radiation therapy is considered to strongly depend on the dosimetric and geometric accuracy of radiation delivery.

Specifically, success is dependent on beam-defining elements such as jaws and MLC leaves accurately following their intended motions, in order to produce the intended modulated fields. Deviations from intended motions can cause dose deviations - doses that are higher or lower than intended in certain locations - and geometric / positioning errors. Dose deviations can create "cold spots" in the target which allow cancer to survive, or excessively "hot spots" that cause tissue necrosis. Geometric errors, which displace radiation from its intended location, can injure adjacent healthy organs. Additionally, because radiation dose is irreversible, patient safety is paramount. Significant excess dose that is delivered to a patient, or dose delivered to the wrong location as a result of a geometric error, cannot be undone and has the potential to severely degrade a patient's quality of life, or be fatal.

For these reasons, a quality assurance program (QAP) is an essential component of radiation therapy. The QAP recognizes that the radiation delivery system, consisting of the accelerator, table, and adjunct imaging and position equipment, is subject to inherent mechanical uncertainties. For example, jaw and MLC leaf positions are accurate only up to certain mechanical tolerances. One goal of the QAP is to quantify the component uncertainties present in the radiation delivery system and, through routine equipment calibration, to maintain those uncertainties at low levels. Residual uncertainties - those which cannot be practically reduced - are accounted for by incorporating safety margins into treatment planning.

For delivery of patient treatment beams, modern radiation therapy attempts to achieve overall geometric accuracy on the order of a few millimeters, and overall dosimetric accuracy on the order of a few percent. A further goal of the QAP is to achieve these end-to-end accuracies. Overall delivery accuracy can be assessed by means of verification imaging.

Verification imaging utilizes a measurement device that produces a quantitative image of a field's dose distribution. The first step in verification imaging is to generate a reference image, i.e., the image that would be measured if delivery were to occur without mechanical deviations or errors. The reference image is often obtained by modeling (predicting) the dose delivered by a patient treatment field to the measurement device. This operation is performed in the treatment planning system or adjunct software. The treatment field is then delivered to the measurement device. The reference and measured images (dose distributions) are compared. Discrepancies between the images that fall outside established tolerance levels represent potential delivery errors, and must be evaluated for their possible impact on the patient dose. Significant discrepancies result in corrective actions, including replanning of the patient's treatment to achieve acceptable agreement between reference and measured doses, and / or repair or calibration of mechanical components, deviations in which contributed to the delivery error.

Verification imaging provides end-to-end verification of the patient treatment beam. Accurate delivery of radiation fields depends on a chain of factors, including correct and accurate commissioning of the treatment planning system (beam modeling); accurate, wi thi -tolerance motion of mechanical components including gantry, jaws, MLC, collimator and table; and accurate dosimetric calibration of the accelerator. Verification imaging verifies that all these components are functioning correctly, to within acceptable tolerance. A general engineering rule states that in order to control a system's accuracy to level X, it is necessary to be able to measure the system to a resolution of XI 10 or better. If measurement resolution is not fine enough, the system's variability will not be recognized, and cannot therefore be controlled. In order to control end-to-end radiation therapy uncertainties to the required levels, it follows that verification imaging must be capable of measuring (detecting) small dose errors (~1 %) and small geometric errors (< lmm).

Verification images can be captured with a variety of two-dimensional measurement devices such as film (Devic, S., "Radiochromic film dosimetry: past, present, and future", Phys. Med. 27, 122-134, 201 1 ) , ion-chamber array (Saminathan, S., Manickam, R., Chandraraj, V, and Supe, S.S., "Dosimetric study of 2D ion chamber array matrix for the modern radiotherapy treatment verification", J. Appl. Clin. Med. Phys. 1 1, 3076, 2010) , diode array

(Gutierrez, A. N., and Calvo, O., "Diode Arrays and QA of Advanced Techniques", Journal of Physic: Conference Series 250, 012049, 2010) , or electronic portal imaging device (EPID) (Van Elmpt, W., McDermott, L.,

Nijsten, S., Wendling, M, Lambin, P., and Mijnheer, B., "A literature review of electronic portal imaging for radiotherapy dosimetry", Radiother. Oncol. 88, 289-309, 2008) . The EPID is a flat-panel array of solid-state electronic sensors, deployed on a movable arm that can be extended from the base of the accelerator gantry, perpendicular to the radiation beam (Fig. 5). Fig. 5 is a side view of a linear accelerator gantry 51, table 52, patient 53, and EPID 54. The EPID 54 extends from the base of the gantry 51 on an articulating arm 55, so that it remains perpendicular to the central axis 56. The EPID 54 can capture a dosimetric image of the delivered beam either wiuh or without the table 52, patient 53 or phantom present. Once extended, beams delivered by the accelerator produce an image of the dose distribution on the EPID 54. One commercial EPID consists of an array of 768 x 1024 pixels, each pixel of size 0.4 mm². Other commercial EPIDs have similar specifications. The spatial resolution of EPIDs is comparable to film, and higher than ion chamber and diode arrays, which have detector spacings of several millimeters. Unlike film, which must be developed and scanned, a further advantage of EPIDs is that their images are captured in electronic form, and automatically stored in a database. Though less common, verification images can also be captured with a three-dimensional measurement device such as a gel dosimeter (De Deene,

Y., "Gel dosimetry for the dose verification of intensity modulated

radiotherapy treatments", Med, Phys. 12, 77-88, 2002) . Although broadly applicable to other measurement devices, including three-dimensional measurement devices, the invention is described here for two-dimension l EPID verification images, examples of which are shown in Fig. 6.

Fig. 6(a) is an example of a reference EPID image. Fig. 6(b) is an example of a measured EPID image containing a delivery error in the circle. The delivery error is too small (i.e., its contrast is too low) to be detectable by eye. The goal of a detection algorithm is to identify such errors, while ignoring deviations caused by wi thin-tolerance motion of beam-defining elements (jaws and MLC leaves). Fig. 6(c) shows the corresponding gamma image (discussed below). Intensity at each point of this image is the computed gamma index. In Fig. 6(d), the corresponding scaled PID image produced by the PID detection algorithm (discussed below). A key component of verification imaging is the detection algorithm that compares reference and measured images, and identifies discrepancies. In Fig. 6, the outlined rectangle in the measured EPID image is a region of the field in which dose is some percentage higher than it should be. The region of higher dose could be caused, e.g., by an MLC leaf or leaves being

significantly out of position during beam delivery. This is an example of a dosimetric error. Alternatively, MLC leaf motions that shape the edge of the field could be incorrect, resulting in a field edge that is shifted from its intended position. This is an example of a geometric error. Geometric and dosimetric errors are closely related, since a geometric error in a field-defining element such as jaw or MLC leaf typically also produces a dosimetric error. The detection algorithm is intended to identify dosimetric and geometric errors caused by out-of-toierance motion of jaws, MLC leaves, gantry, collimator, etc.

Desirable properties of the detection process are as follows:

1) The process should exhibit high sensitivity. It must be capable of detecting small (~1%, <lmm) dosimetric and geometric errors caused by out- of-tolerance motion of jaws and MLC leaves.

2) The process should exhibit high specificity. A poorly designed process will incorrectly interpret image noise and residual uncertainties as delivery errors. This has the effect of producing large numbers of false positives— supposed delivery errors that consume time and effort for investigation, but turn out to be caused by within- tolerance motion of jaws, MLC leaves, etc. The process should ignore discrepancies caused by within-tolerance motion, but correctly identify errors caused by out-of-tolerance motion.

3) The process should incorporate a method to identify levels of noise and residual uncertainties in verification images, so that it can be calibrated to ignore within-tolerance motions of jaws, MLC leaves, etc.

4) The process should incorporate a method to quantify detection performance in terms of the standard statistical metrics of sensitivity, specificity and overall accuracy at detecting delivery errors. Practically, such a method allows detection algorithm parameters to be calibrated to produce a specified level of detection accuracy. In general, it is undesirable to use a detection process whose detection accuracy has not been quantified.

The most widely used detection algorithm is the gamma algorithm proposed by Low et al (Low, D. A., Harms, W. B., Mutic, S., and Purdy, J. A., "A technique for the quantitative evaluation of dose distributions", Med. Phys. 25, 656-661 , 1998) . The gamma algorithm associates a real -valued gamma index γ with each pixel or voxel u in the measured image. The gamma index is given by the formula (Low et al,, 1998)

where Ι,„, is the image intensity (dose) at pixel (voxel) u in the measured image, is the image intensity at pixel (voxel) v in the reference image, d(u,v) is the distance between u and v (2D or 3D as appropriate), f_d0Se is the γ dose criterion, f₍n_sl is the γ distance criterion, and the minimum is over a set of reference image positions V, typically the set of pixels (voxels) v falling within some radius of u. Low et al proposed that the nomialization intensity ionn be set equal to the image's maximum intensity, although it can also be set equal to the local intensity (i.e., intensity at u or v). The derived image consisting of gamma index values is referred to as the gamma image. An example of EPID verification images and the corresponding gamma image is given in Fig. 6.

For every pixel u in the measured image, the intention of the gamma algorithm is to find the closest neighboring pixel v in the reference image that best matches the dose at u. The parameters fd_0se and f ist are the gamma algorithm's intended tolerance to dosimetric and geometric errors. The first term under the square root in equation (1) is the percent dose difference divided y fdose · The second term is the positional offset divided by fd_tst · Combinations of dose differences greater than f _ose and positional offsets greater than f^, produce gamma values greater than one. Conversely, if the dose difference and positional offset are both small with respect to fdose and fdist , the gamma value will be less than one. Image classification is based on the percent of gamma indices below a threshold of one. The measured image is classified as error- free if the percent of gamma values below one exceeds a cutoff τ. This criterion is illustrated in Fig. 7 which is an illustration of image classification according to the gamma algorithm. Gamma indices are computed for all pixels and their probability distribution is computed. Image classification is based on the percent of gamma values less than one. If this percentage exceeds a cutoff (e.g., 95%) the image is classified as error-free. Gamma indices above one are indicative of a delivery error.

Undesirable aspects of the gamma algorithm are as follows:

(a) After it was proposed, the algorithm was recognized as being sensitive to the statistical noise present in images. (Low, D. A., and Dempsey, J, F., "Evaluation of the gamma dose distribution comparison method", Med. Phys. 30, 2455-2464, 2003). Depending on how the gamma algorithm is implemented, it can be shown that the noise present in EPID images reduces the sensitivity, specificity and detection accuracy of the algoritlim. (Gordon, J. J., Gardner, J. K., Wang. S., and Siebers, J. V., "Reliable detection of fluence anomalies in EPID-based IMRT pretreatment quality assurance using pixel intensity deviations", Med. Phys. 39, 4959-4975, 2012). It is possible this observation applies also to other high resolution detectors, such as film.

Commercial implementations of the gamma algorithm do not explicitly account for noise.

(b) The intention of the gamma algorithm is to detect dosimetric or geometric errors greater t &n fd_ose or fdist, ^an<l ignore those less than i_0Se or fdist . However, the algorithm provides no methodology for verifying detection performance. For EPID images, rigorous analysis shows that the gamma algorithm can fail to detect errors greater than fd_0Se or f i_st, due to interaction with image noise and / or poor implementation. (Gordon et al., 2012) It is desirable that there be some method of quantifying the algorithm's detection performance, in terms of sensitivity, specificity and accuracy,

(c) The gamma algorithm does not include any methodology for setting the values of parameters fdose , fdist and τ, which control the algorithm's sensitivity, specificity and detection accuracy. Users typically choose values that are commonly used by others, e.g., f_dose = 3%, fdi_sS ~ 3mm and τ = 95%.

Logically, detection performance will depend on the properties of the verification images, including the level of noise, and the level of residual positional uncertainties induced by mechanical deviations of field-defining elements. It is desirable to have a prescription for setting algorithm parameters based on measured image properties.

The search for a minimum in equation (1)— referred to as a distance-to-agreement (DTA) search— is numerically expensive.

Consequently, several variants of the gamma algorithm have been proposed that produce similar gamma indices, but with less computation. Suppose one maps n-dimensional verification images to an (n+1 )-dimensional space, where the (n+l)st dimension is image intensity. (In the case of 2-dimensional EPID images, x and y dimensions would be pixel position, and the z-dimension would be image intensity.) The parameters fdose and fdist define an ellipsoid region around each pixel in the (n+l)-dimensional space, within which the gamma value is less than one. As long as the reference image dose surface intersects the ellipsoid, the minimum gamma value produced by the DTA search will be less than one. This is illustrated for a one-dimensional image in Fig. 8a which is an illustration of one-dimensional profiles through reference and measured verification images. Image intensity (dose) is plotted on the y- axis against pixel number or distance on the x-axis. In the gamma algorithm, parameters f_dose and fdhi define an ellipsoid region around each pixel, within which the gamma index is one. Taken over all pixels in the measured image, the envelope of ellipsoids produces an acceptance "tube", illustrated in Fig. 8b. Fig. 8b shows the surface created by the envelope of ellipsoids is an acceptance tube. Provided the reference image dose surface remains wholly within the acceptance tube, the measured image is classified as error-free.

Bakai et al (Bakai, A., Alber, M., and Nusslin, F., "A revision of the gamma-evaluation concept for the comparison of dose distributions", Phys. Med. Biol, 48, 3542-3553, 2003) and Blanpain and Mercier (Blanpain, B., and

Mercier, D., "The delta envelope: a technique for dose distribution

comparison", Med. Phys. 36, 797-808, 2009) proposed variants of the gamma algorithm, which are based on direct computation of the ellipsoid acceptance tube. Computation of this tube exploits knowledge of the local dose gradient at each pixel. These algorithms can therefore be considered to scale the gamma index at a point in the measured image by the local gradient at that point. The same idea underlies the approach of Moran et al (Moran, J. M. Radawski, J., and Fraass, B. A., "A dose gradient analysis tool for IMRT QA", Journal of Applied Clinical Medical Physics 6, 2005), which directly scales the dose difference at a point by the dose gradient at that point.

Importantly, none of these approaches scales the gamma index by a factor that depends on image noise and the positional variability of field-defining elements, which is a feature of the present invention. They perform a simpler scaling based on dose gradient only. The gamma variants noted above suffer from the same disadvantages as the original gamma algorithm: they do not account for image noise; they do not provide a method of setting detection parameters based on measured image properties, such as residual positional uncertainties; and they do not provide a method of quantifying the achieved detection performance. SUMMARY OF THE INVENTION

The invention comprises a process for detecting delivery errors in a measured image of a radiation therapy treatment field. A reference image is independently obtained, assuming error-free delivery. For each pixel or voxel, the percent intensity difference between measured and reference images is computed. Each percent intensity difference is scaled by a quantity that depends on the local intensity gradient, the measured image noise, and the measured positional variability of field-defining elements. Optionally, further processing such as median filtering is applied to the derived image consisting of scaled percent intensity differences. Image classification is based on the percent of resulting scaled pixel values below a threshold. The measured image is classified as error-free if the percent of scaled pixel values below the threshold exceeds a cutoff

Noise and the variability of field-defining elements are measured either in a single image, or in a population of images, by plotting the variance of percent intensity differences against the directional intensity gradients squared, and measuring the intercept and slopes of these relationships. The invention further comprises a method for quantifying detection accuracy. This method inserts simulated errors into measured images and quantifies sensitivity, specificity and detection accuracy. This method is used to identify values for the threshold and cutoff that produce acceptable detection accuracy. Novel elements of the invention are the method of scaling percent intensity differences, the method of measuring image noise and the positional variability of field-defining elements, and the method of quantifying and calibrating detection performance. BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing and other objects, aspects and advantages will be better understood from the following detailed description of a preferred embodiment of the invention with reference to the drawings, in which:

Fig. 1 is illustration of a linear accelerator, showing table and gantry rotation;

Fig. 2 is a section through the linear accelerator central axis, showing beam defining elements: jaws and multileaf collimator (MLC) leaves;

Fig. 3 is a beams-eye view from the linear accelerator source along the central axis showing an irregular field shaped by jaws and MLC leaves;

Fig. 4 is an example of a modulated treatment field, red indicating regions of high radiation dose and blue indicating regions of low dose;

Fig. 5 is a side view of a linear accelerator, showing an electronic portal dosimetry device extended from the base of the gantry so as to capture a dosimetric image of the delivered field;

Fig. 6 shows examples of a reference EPID image (a), the

corresponding measured image (b) showing an area in which a dosimetric error has occurred (circled), the corresponding gamma image (c), and the corresponding PID image (d);

Fig. 7 shows an example of the gamma index distribution from a measured EPID image;

Fig. 8a shows the acceptance ellipsoid defined by parameters fdose and fdist for the gamma algorithm;

Fig. 8b shows the acceptance "tube" constructed as the envelope of acceptance ellipsoids;

Figures 9a and 9b, taken together, are a flow diagram illustrating the detection process according to the invention;

Fig. 10 illustrates image classification using the PID process according to the invention;

Figures 1 la and 1 lb, taken together, are a flow diagram illustrating the method of quantifying noise and positional deviations according to the invention;

Fig. 12 illustrates plots of PID variance versus gradient squared, used to determine variances of positional deviation of field defining elements; and

Figures 13a and 13b, taken together, are a flow diagram illustrating the method of quantifying detection performance and optimization of κ and τ according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

A concrete implementation of the invention is described for two- dimensional EPID images. The invention generalizes to other two- dimensional images, including those produced by other measurement devices, and to three-dimensional images. In the case of three-dimensional images, standard deviations of positional deviations must be estimated for three instead of two dimensions, and the scaling factor for percent dose differences should include a term for the third dimension. The following assumes that the EPID has been calibrated according to site policy or manufacturer

instructions.

The invention comprises: (i) a detection process; (ii) a method for quantifying noise and standard deviations of positional deviations, which are required as parameters in the detection process, and (iii) a method for quantifying detection performance, which allows thresholds and cutoffs employed in the detection process to be optimized. The invention is implemented in computer software. Detection Process Reference is now made to Figs. 9a and 9b which illustrate the detection process according to the invention. For a treatment field, a reference EPID image is generated at step 901. The method by which this is done is well known in the prior art and therefor outside the scope of the invention. The resulting image consists of an N_x x N_y array of scalar intensities, where N_x is the number of pixels in the x (crossplane) direction and N_y is the number of pixels in the y (inplane) direction. Fig. 6a shows an example.

At step 902, at the linear accelerator, the treatment field is delivered using the same EPID acquisition mode and under the same conditions as were assumed for the reference image, and a measured EPID image is acquired.

Optionally, the measured image is multiplied by a scale factor at step

903 to compensate for linear accelerator output variations. This can be done according to a variety of methods. One method is to identify pixels in the reference image having intensity greater than e.g., 20% of maximum intensity, then scale the measured image so that the average intensity over this subset of pixels agrees with the reference image. The measured image is classified as errored if the scale factor differs from one by more than a specified tolerance that is established through statistical analysis of accumulated data, or arbitrarily set (e.g., to 2%).

Optionally, at step 904, the measured image is aligned to the reference image using translation, rotation and uniform scaling. The goal is to obtain best alignment between the fields in each image. Edges of the measured image are trimmed and / or padded to the original size. The measured image is classified as errored if corrections differ from baseline values by more than a specified tolerance, that are established through statistical analysis of accumulated data, or arbitrarily set. The baseline values for translational and rotational corrections are expected to be close to zero, non-zero values reflecting effects such as gantry and EPID sag at non-zero gantry angles. The baseline value for scale corrections is expected to be close to one, non-one values reflecting small deviations in the source-to-imager distance.

Optionally, image regions encompassing the treatment beam are identified in step 905 in the reference and measured images, for example, by identifying pixels with intensity greater than or equal to 50% of maximum intensity. An overlap coefficient, such as the Dice overlap coefficient, is computed for the two regions. The measured image is classified as errored if the overlap coefficient differs from a baseline value (e.g., one) by more than a specified tolerance, that is established through statistical analysis of accumulated data, or arbitrarily set.

The percent intensity difference (PID) is calculated at step 906 for each pixel in the measured image, according to the formula: δ,· = 100 ^■ (p_j—p_j ) l P_j , where p . and p_} are theyfh pixel intensities in the reference and measured images. A PID image is generated, such that the intensity at each pixel is the corresponding PID.

A scaled PID image is generated at step 907, by scaling each PID value according to the formula:

where R_XJ- is the relative image gradient (percent intensity change per pixel) in the reference image's x-direction at pixel j, R_yj is the relative image gradient in the reference image's y-direction at pixel j, σ# is the standard deviation of image noise, σ_Λ is the standard deviation of position deviations of beam defining elements in the x-direction, and _y is the standard deviation of position deviations of beam defining elements in the y-direction. R_xj and R_v are calculated from the reference image, σ_# , σ_ν and a_y are calculated according to the accompanying "Method of quantifying noise and positional deviations". An example of a scaled PID image is shown in Fig. 6d.

Optional filtering is applied to the scaled PID image at step 908. This could utilize a variety of established image processing techniques, designed to emphasize image features in the presence of noise. One example would be median filtering.

A subset of pixels is identified at step 909, representing those pixels of practical interest for delivery error detection. This could be done in a variety of ways. A suggested approach would be to identify all pixels that have intensity greater than equal to 50% of maximum dose in the reference image and include all pixels that are within this set, or within some maximum distance (e.g., 5mm) from a pixel in this set. The intent is to exclude the majority of pixels in the periphery of the image, that are distant from the treatment field and so irrelevant to error detection.

The distribution of processed PID values is computed at step 910 for the subset of pixels identified in step 909. The fraction of these pixels with a processed PID value that is less than κ standard deviations from the mean value is computed: φ = PR[ |δ- μ| <κσ]. At step 91 1, the measured image is classified as error-free if φ > τ and errored if φ < τ, where τ is a cutoff. The classification process is illustrated in Fig. 10. Here, image classification is performed according to the PID process of the invention. Scaled (and optionally processed) PID values for most pixels are expected to fall in the range [-κσ,+κσ]. Due to statistical variability, a small number of pixels may have PID values outside this range, even though there is no delivery error.

Significant numbers of pixels having scaled PID values outside [-κσ,+ σ] are indicative of a delivery error.

Values for parameters κ and τ can be predefined. A reasonable approach would be to set κ = 5, and set τ to reflect the minimum size of error that one is attempting to detect. For example, if one is attempting to detect delivery errors covering regions greater than 10^x 10 = 100 pixels, and the subset of pixels identified in step 909 contains 10,000 pixels, one would set τ = 1 - (100/10000) = 0.99. Values for κ and τ that may give better detection performance for a range of error sizes can be identified using the

accompanying "Method of quantifying detection performance and

optimization of K and τ".

Method of quantifying noise and positional deviations

Referring now to Figs. 1 l a and 1 lb, there is shown the method of quantifying noise and positional deviations. At step 1 101 , analysis may be performed on a single treatment field (reference image) for which one or more repeated measured images have been acquired, or multiple treatment fields (reference images) for which one or more repeated measured images have been acquired. Optionally, at step 1 102, scale factors are applied to measured images as in step 903 of Fig. 9a.. Optionally reject (i.e., exclude from the analysis) those images for which the scale factor differs from one by more than an acceptable tolerance.

Optionally, at step 1 103, measured and reference images are aligned as in step 904 of Fig. 9a.. Optionally reject those images for which corrections differ from baseline values by more than an acceptable tolerance. Optionallly, at step 1 104, field regions are identified as in step 905 in Fig. 9a. Optionally reject those images for which overlap coefficients differ from baseline values by more than an acceptable tolerance.

At step 1105, PIDs are calculated as in step 906 in Fig. 9b, The following steps are applied to the combined set of PIDs from all surviving measured images.

In the reference image, the relative gradients Rxj and Ry j in the x- and y-directions, and the total gradient R = -^R + ^_y- , are computed for each pixel at step 1 106.

At step 1 107, σ_# is computed as the standard deviation of PID values for pixels where R_/ is close to zero, e.g., Ry < 1.

PID values are allocated to bins at step 1 108 based on the x-gradient R² _{X }} . The total range of squared gradient values are divided into a number of bins. At step 1 109, the variance of the PIDs in each bin is plotted against i?_t ² where R² is the center value of the bin. An example of this plot is the plot labeled "OVCR R_x" in Fig. 12. A straight line is fitted to the plot and detection algorithm parameter σ_χ is set equal to the square root of the line's slope.

At step 1 110, PID values are allocated to bins based on the y-gradient R_y squared. The total range of squared gradient values are divided into a number of bins. At step 1 11 1, the variance of the PIDs in each bin is plotted against R² , where R² is the center value of the bin. An example of this plot is the plot labeled "OVCR R_y" in Fig. 12. A straight line is fitted to the plot and detection algorithm parameter σ_γ is set equal to the square root of the line's slope. Fig. 12 shows plots of the PID variance versus relative x-gradient squared (labeled "OVER R_x") and relative y-gradient squared (labeled "OVER R_y"). Data was generated by binning pixels according to x- or y- gradient, and computing the variance of the corresponding PID values. The gray plots are fitted lines.

Method of quantifying detection performance and optimization of κ and τ

Reference is next made to Figs. 13a and 13b which illustrate the method of quantifying detection performance and optimization of and τ. At step 1301 , analysis may be performed on a single treatment field (reference image) for which one or more repeated measured images have been acquired, or multiple treatment fields (reference images) for which one or more repeated measured images have been acquired. Optionally, at step 1302, scale factors are applied to measured images as in step 903 in Fig. 9a. Optionally, those images for which the scale factor differs from one by more than an acceptable tolerance are rejected (i.e., excluded from the analysis).

Optionally, measured and reference images are aligned at step 1303 as in step 904 in Fig. 9a. Optionally, those images for which corrections differ from baseline values by more than an acceptable tolerance are rejected.

Optionally, field regions are identified at step 1304 as in step 905 in Fig. 9a. Optionally, those images for which overlap coefficients differ from baseline values by more than an acceptable tolerance are rejected. Randomly, at step 1305, select one of the original measured images, i.e., prior to steps 1302, 1303, 1304, corresponding to an image in the surviving pool. The selected image is classified as error-free or errored using the detection process of Figs. 9a and 9b. For each image that is classified as error-free, increment the true negative count. For each image that is classified as errored, increment the false positive count. At step 1306, step 1305 is repeated N_neg times. True negative and false positive counts are divided by N_neg to obtain true negative and false positive rates.

At step 1307, one of the original measured images is randomly selected. A subset of pixels that are of practical interest for error detection is identified, as in step 909 in Fig. 9b. A simulated delivery error is randomly inserted at a random location in this set of pixels, in such as way that errored pixels fall wholly within the set of identified pixels. The simulated delivery error can take any form. For example, it could consist of a region of n x m pixels in which intensity (dose) is raised or lowered by q%, where parameters n, m and q are predefined. The selected image is classified as error- free or errored using the detection process (Figs. 9a and 9b). For each image that is classified as error-free, increment the true positive count. For each image that is classified as errored, increment the false negative count.

At step 1308, step 1307 is repeated N_pos times. True positive and false negative counts are divided by N_pcis to obtain true positive and false negative rates. The sum of true negative and true positive counts is divided by (N_neg + N_p0S ) to obtain detection accuracy.

For a given delivery error, compute in step 1309 detection accuracy (steps 1301 - 1308) for a fixed value of the threshold parameter κ and a range of values of the cutoff parameter τ. Identify the value of τ that provides maximum detection accuracy. Incrementally vary κ and repeat the foregoing procedure to find the combination of κ and τ that provides maximum detection accuracy. At step 1310, step 1309 is repeated for the range of simulated delivery errors that one wishes to detect. Find the combinations of and τ that provide maximum detection accuracy for specific simulated errors.

While the invention has been described in terms of a single preferred embodiment, those skilled in the art will recognize that the invention can be practiced with modification within the spirit and scope of the appended claims.

Claims

1. A process for detecting delivery errors in a measured image of a radiation therapy treatment field, comprising the steps of:

independently obtaining a reference image using a first mode of an imaging device;

delivering a plurality of treatment fields using the first mode and under the same conditions as were assumed for obtaining the reference image, and acquiring a measured image for each treatment field;

computing a percent intensity difference (PID) between measured and reference images;

scaling the percent intensity difference by a value comprising a local dose gradient, measured image noise, and measured positional variability of field-defining elements; and

classifying images based on the scaled percent intensity difference based on a threshold, wherein the measured image is classified as error-free.

2. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 1 , further comprising the step, prior to PID computation, of multiplying the measured image by a scale factor to compensate for linear accelerator output variations and classifying the measured image as errored if the scale factor differs from one by more than a specified tolerance.

3. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 1 , further comprising the step, prior to PID computation, of aligning the measured image to the reference image using translation, rotation and uniform scaling to obtain best alignment between the fields in each image and classifying the measured image as errored if corrections differ from baseline values by more than a specified tolerance.

4. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 1, further comprising the step, prior to PID computation, of identifying image regions encompassing the treatment beam in the reference and measured images, and classifying the measured image as errored if an overlap coefficient differs from a baseline value by more than a specified tolerance.

5. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 1, further comprising the step, after the computation of scaled PIDs, of applying image processing techniques to the measured image consisting of scaled percent intensity differences to emphasize delivery errors.

6.The process for detecting delivery errors in a measured dosimetric image of a radiation therapy treatment field of claim 1 , which process comprises a detection process comprising the steps of:

obtaining a reference image, the image consisting of a two- dimensional or three-dimensional array of scalar values, the scalar values being image intensities;

delivering a treatment field using the same acquisition mode and under the same conditions as were assumed for the reference image, and acquiring a measured image;

multiplying the measured image by a scale factor to compensate for linear accelerator output variations and classifying the measured image as errored if the scale factor differs from one by more than a specified tolerance; aligning the measured image to the reference image using translation, rotation and uniform scaling to obtain best alignment between the fields in each image and classifying the measured image as en red if corrections differ from baseline values by more than a specified tolerance;

identifying image regions encompassing the treatment beam in the reference and measured images, and classifying the measured image as errored if an overlap coefficient differs from a baseline value by more than a specified tolerance;

calculating percent intensity difference (PID) for each pixel in the measured image, according to the formula: 5j = 100 (p_j - p_} ) I p . , where p_} and p_j are the jt pixel intensities in the reference and measured images and generating a PID image, such that the intensity at each pixel is the

corresponding PID;

generating a scaled PID image by scaling each PID value according to the formula:

where R_xj is the relative image gradient (percent change per pixel) in the reference image's x-direction at pixel j, R_yj is the relative image gradient in the reference image's y-direction at pixel j, ϋβ is the standard deviation of image noise, σ_χ is the standard deviation of position deviations of beam defining elements in the x-direction, and a_y is the standard deviation of position deviations of beam defining elements in the y-direction, R_x and R_Vi are calculated from the reference image;

identifying a subset of pixels representing those pixels of practical interest for delivery error detection to exclude a majority of pixels in a periphery of the image that are distant from the treatment field and so irrelevant to error detection;

computing a distribution of processed PID values for the subset of pixels; and computing a fraction of these pixels with a processed PID value that is less than κ standard deviations from the mean value as φ = PR[ |δ-μ| <κσ ] and classifying the measured image as error-free if φ > τ and errored if φ < τ, where τ is a cutoff.

7. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 1. which process comprises a method of quantifying noise and positional deviations comprising the steps of: analyzing one or more reference images for which one or more repeated measured images have been acquired;

calculating percent intensity difference (PID) for each pixel in each measured image;

computing in the reference images, relative gradients R_xj and R_yj in the x- and y-directions, and the total gradient R_j = jR _j + R , , for each pixel; computing o_B as the standard deviation of PID values for pixels where

R_j is close to zero, e.g., R_j < 1 ;

allocating PID values to bins based on the x-gradient . and dividing a total range of squared gradient values into a predetermined number of bins;

plotting a variance of the PIDs in each bin against

, where R_x ² is the center value of the bin, fitting a straight line to the plot and setting σ_χ equal to the square root of the line's slope; and

allocating PID values to bins based on the y-gradient R_yJ squared and dividing a total range of squared gradient values into a predetermined number of bins; and

plotting a variance of the PIDs in each bin against R_y ² , where Ic_y is the center value of the bin, fitting a straight line to the plot and setting σΥ equal to the square root of the line's slope.

8. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 7, wherein the process for quantifying noise and positional deviations further comprises the step, prior to

PID computation, of multiplying each measured image by a scale factor to compensate for linear accelerator output variations and rejecting by omitting from further analysis the measured image if the scale factor differs from one by more than a specified tolerance.

9. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 7, wherein the process for quantifying noise and positional deviations further comprises the step, prior to PID computation, of aligning the measured image to the reference image using translation, rotation and uniform scaling to obtain best alignment between the fields in each image and rejecting by omitting from further analysis the measured image if the scale factor differs from one by more than a specified tolerance.

10. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 7, wherein the process for quantifying noise and positional deviations further comprises the step, prior to PID computation, of identifying image regions encompassing a treatment beam in the reference and measured images and rejecting by omitting from further analysis the measured image if an overlap coefficient differs from a baseline value by more than a specified tolerance.

1 1. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 7, wherein the process for quantifying noise and positional deviations further comprises the steps of, prior to PID computation, of identifying a subset of pixels representing those pixels of practical interest for delivery error detection to exclude a majority of pixels in a periphery of the image that are distant from the treatment field and so irrelevant to error detection, and obtaining gradient plots only for the identified subset of pixel.

12. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 1 , which process comprises a method of quantifying detection performance comprising the steps of:

analyzing one or more reference images for which one or more measured images have been acquired;

for Nncg trials, randomly selecting one of the original measured images and classifying the selected image as error-free or errored using the detection process;

for each image that is classified as error-free, incrementing a true negative count, and for each image that is classified as errored, incrementing a false positive count;

dividing negative and false positive counts by N„_eg to obtain true negative and false positive rates;

for Npos trials, randomly selecting one of the original measured images, identifying a subset of pixels that are of practical interest for error detection, randomly inserting a simulated delivery error at a randomly selected location in this set of pixels, and classifying the selected image as error-free or errored using the detection process;

for each image that is classified as errored, incremening the true positive count, and for each image that is classified as error-free, incrementing the false negative count; and

dividing true positive and false negative counts by N_pos to obtain true positive and false negative rates, and dividing a sum of true negative and true positive counts by (N_neg + N_pos ) to obtain detection accuracy.

13. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 12, wherein the method for quantifying detection performance further comprises the step of multiplying each measured image by a scale factor to compensate for linear accelerator output variations and rejecting by omitting from further analysis the measured image if the scale factor differs from one by more than a specified tolerance.

14. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 12, wherein the method for quantifying detection performance further comprises the step of aligning the measured image to the reference image using translation, rotation and uniform scaling to obtain best alignment between the fields in each image and rejecting by omitting from further analysis the measured image if the scale factor differs from one by more than a specified tolerance.

15. The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 12, wherein the method for quantifying detection performance further comprises the step of identifying image regions encompassing a treatment beam in the reference and measured images and rejecting by omitting from further analysis the measured image if an overlap coefficient differs from a baseline value by more than a specified tolerance.

16 The process for detecting delivery errors in a measured image of a radiation therapy treatment field of claim 12, which method comprises a method of optimizing detection performance comprising the steps of: for a given delivery error, computing detection accuracy for a fixed value of the threshold parameter and a range of values of the cutoff parameter τ to identify a value of τ that provides maximum detection accuracy; and

incrementally varying κ and repeat the foregoing procedure to find the combination of κ and τ that provides maximum detection accuracy for specific simulated errors.

17. A method and statistical validation technique for detecting differences between radiation therapy verification images in order to detect delivery errors in radiation therapy treatment fields, comprising the steps of:

utilizing verification images in which intensity is proportional to dose based on comparison of reference and measured images;

detecting and correcting output variations between the images, classifying the measured image as errored if the output variation exceeds a predetermined value;

aligning the images using translations, rotations and scale changes, and classifying the measured image as errored if the corrections exceed a predetermined value;

comparing field outlines in the measured and reference images by means of an overlap coefficient, and classifying the measured image as errored if the overlap is below a predetermined value;

computing a percent intensity difference (PID) at each pixel;

scaling each PID by a factor according to the formula:

where R_XJ- is the relative image gradient (percent change per pixel) in the reference image's x-direction at pixel j, R_yj is the relative image gradient in the reference image's y-direction at pixel j, σ# is the standard deviation of image noise, σ_χ is the standard deviation of position deviations of beam defining elements in the x-direction, and a_y is the standard deviation of position deviations of beam defining elements in the y-direction. R_xj and R_yj- are calculated from the reference image;

presenting an array of scaled PIDs in the form of an image;

identifying a subset of pixels which are to be included in error detection;

computing a distribution of scaled PEDs for the identified subset of pixels and determining a fraction φ lying less than κ standard deviations from the mean, where κ is a parameter; and

classifying the measured image as error- free if op > τ or errored if φ<τ , where τ is a parameter.

18. A method of quantifying detection performance for image analysis comprising the steps of:

for N_neg trials, randomly selecting one of the original measured images and classifying the selected image as error-free or errored using the detection process;

dividing negative and false positive counts by N_neg to obtain true negative and false positive rates;