fig1 shows a radiotherapy system 1 according to embodiments of the present invention . structurally , the system is similar to a conventional radiotherapy apparatus . the system comprises a gantry 10 on which is mounted a source of radiation 12 and , diametrically opposite the source 12 , a radiation detector 16 . such detectors are commonly referred to as portal imagers . the radiation source 12 is typically a linear accelerator producing x - rays or other penetrating radiation . the gantry is rotatable about an axis 22 . in the figures , the gantry 10 is depicted as a ring - shaped support . alternatively , however , the gantry may comprise a c - arm , with the source 12 and imager 16 on opposite arms . the isocentre of the system is defined as a plane running through the rotation axis 22 of the gantry 10 perpendicular to the instantaneous axis of the radiation beam . a collimator 14 is coupled to the radiation source 12 in order to collimate and shape the radiation beam . that is , a first collimation of the radiation ( not illustrated ) takes place close to the source 12 . this collimates the radiation produced by the source into a beam , e . g . a cone - or fan - shaped beam diverging away from the source . a further collimator 14 then acts on this collimated beam in order to shape the radiation as required for therapy . an example of a suitable collimator for this aspect is a multi - leaf collimator ( mlc ). such devices comprise one or more banks of parallel leaves , each of which can be moved in a direction transverse to the radiation beam axis . the leaves are moveable into and out of the path of the radiation beam to a greater or lesser extent , and thus the combination of leaf positions collectively defines a shaped aperture through which radiation passes . in one embodiment , the mlc comprises two banks of leaves positioned on opposite sides of the radiation beam , with each leaf able to take any position with a range from outside the radiation beam to crossing the radiation beam . in order to sufficiently attenuate ( i . e . block ) the high - energy radiation , the leaves have a significant depth in a direction parallel to the radiation beam axis , and are manufactured from high atomic number materials such as tungsten . thus , the output of the radiation source 12 and the collimator 14 is a shaped radiation beam 24 directed generally inwards towards the axis of rotation of the gantry . control and processing circuitry 26 is in communication with the gantry 10 , the source 12 and the collimator 14 and controls their operation . a patient 20 is positioned on a support 18 for treatment such that a treatment target 21 ( e . g . a tumour ) is placed at the isocentre of the system . the longitudinal axis of the support and , thus , the patient 20 usually but not necessarily lie substantially parallel to the rotation axis 22 of the gantry . various processes and apparatus for positioning and locating the patient will be familiar to those skilled in the art . in one embodiment , the support allows linear translation of the patient in three dimensions . the support 18 may also allow for tilting and rotation of the patient 20 , thus providing movement in up to six degrees of freedom ( i . e . x , y , z , pitch , yaw and roll ). in operation , the system 1 performs the method as set out in fig2 . the method begins once the patient has been positioned correctly , i . e . with the target at the isocentre . embodiments of the invention may also provide for motion of the patient during treatment . for example , the support 18 may compensate for movement of the patient due to the respiratory or cardiac cycles ( i . e . minimizing the motion of the target 21 relative to the system 1 ). similarly , the movement and positioning of the collimator leaves may compensate for such cyclical movement of the target 21 . in the illustrated method , however , such processes are not considered for simplicity . in step 100 , the method begins with the radiation source 12 generating a radiation beam . the radiation itself may be x - rays or other penetrating radiation as required . in step 102 , the radiation beam is collimated by action of the collimator 14 to generate a shaped , collimated beam 24 for the purposes of therapy . for example , the beam may be shaped to conform to the shape of the target , or a part of the target , or any other shape to achieve a desired dose distribution in and around the target . this beam is incident on the patient 20 and , inevitably , some of the radiation is absorbed . the amount of absorption depends on the particular structures within the patient . in step 104 , the attenuated radiation beam is detected by the portal imager 16 . the detection signals are converted to data and provided to the control circuitry 26 . these data sets largely comprise pixel data and can therefore ( usually ) be assembled into viewable images . they will therefore be referred to as “ images ” in this application , but this should not be interpreted to mean that the data sets must be stored in an image format , or that they must be viewed in an image form at any point . often , the data sets will be transferred to the processor for conversion into ct datasets without ever being viewed in an image form . steps 100 , 102 and 104 generally all take place while the gantry 10 rotates around the patient ( step 106 ). thus steps 100 , 102 , 104 and 106 may all take place substantially simultaneously . that is , the gantry 10 rotates around the patient , while the source 12 continuously generates a radiation beam and the collimator leaves move to new positions . the control circuitry 26 controls this operation , and may vary the radiation energy , the dose rate , the collimator positions and the speed of rotation of the gantry throughout the process . the method then proceeds to step 108 in which the circuitry 26 applies an algorithm to the portal image data and reconstructs a three - dimensional ct image of the target region . there is no fixed time at which the algorithm is applied . the algorithm acts on the data that has been acquired up until the point the algorithm is invoked . if a three - dimensional ct image is required during treatment , the algorithm may be invoked during treatment . the ct image can then be used to guide the therapy for the remainder of the treatment . alternatively , the ct image may be generated after treatment to assist the determination of the dosage distribution delivered to the patient , or to record the internal motion of the patient &# 39 ; s tissues . generally , however , the more data that is acquired , the better the quality of the ct image that is obtained . the ct image so generated is a measure of the attenuation at any one point in three - dimensions within the target region . thus , in one embodiment , the attenuation coefficient itself ƒ ( r , φ , z ) is found ( where ( r , φ , z ) are cylindrical polar coordinates ). in alternative embodiments , the variable that is reconstructed may be an alternative quantity that is nonetheless related to the attenuation coefficient and provides useful imaging information . hereinafter , the term “ attenuation coefficient ” is taken to mean the attenuation coefficient ƒ ( r , φ , z ) as well as these related quantities . the algorithm that is applied in step 108 allows the generation of a three - dimensional image from portal images despite many challenges . the portal images have a very narrow field of view of the target , as the collimator leaves act to block any radiation that is not directed towards the target in order to minimize the damage to surrounding healthy tissue . moreover , the positions of the collimator leaves may change as the gantry rotates around the patient . thus , the field of view at one angle of rotation will not in general be the same as the field of view in another angle of rotation . the algorithm has several mathematical components , each of which act — both alone and in combination — to overcome these difficulties . a first part of the algorithm relates to an acknowledgement that it may not be possible to reconstruct the attenuation coefficient for all locations within the target region . first , it is impossible to image locations that fall outside the target region entirely — and hence the path of the radiation beam at all gantry rotation angles . there simply is no data to reconstruct values from . however , there may be other locations within the target region which only fall within the path of the radiation beam for a subset of the gantry rotation angles . for these locations , there may nonetheless be insufficient data to reconstruct an accurate value for the attenuation coefficient . reconstructing values of the attenuation coefficient for these locations may reduce the overall quality of the image . thus , for each location within the target region , the first part of the algorithm sets the attenuation coefficient to a null value if that location falls outside the radiation beam for more than a threshold angular extent . for example , that threshold may be set at 270 °. in that case , if a location lies inside the radiation beam for only 75 ° of the 360 ° revolution of the gantry , the attenuation coefficient for that location is set to null . the null value may be zero or any other value which in practice is regarded by the system as null . a second part of the algorithm normalizes the reconstructed value of the attenuation coefficient for a particular location according to the angular range that location was inside the radiation beam . thus , if a particular location lies within the radiation beam for a fraction of a complete revolution , the reconstructed attenuation coefficient for that location is divided by that fraction . for example , if the location lies within the radiation beam for 180 ° ( i . e . half of the possible gantry angles ), the value of the attenuation coefficient may be multiplied by two to account for this . both of these parts of the algorithm may employ a masking function to define the position of the collimator leaves at any particular rotation angle . in a plane orthogonal to the radiation beam axis , the masking function is equal to zero for locations that are blocked by the collimator leaves and one for locations that are open for radiation to travel through . the two variables used to define the location in that plane may be parallel and orthogonal respectively to the direction of travel of the collimator leaves , or at any angle with respect to the leaves . part of the algorithm is a measure of the energy received at the portal imager 16 . previously we have described that it may not be possible to reconstruct values of the attenuation coefficient for locations that lie outside , or partially outside , the radiation beam during rotation of the gantry . however , in order to reconstruct values for locations that lie inside the radiation beam , it may nonetheless be necessary to estimate values for the energy that would have been received by the portal imager for those locations that lie outside the radiation beam . that is , “ would have been received ” had the collimator leaves , for example , not blocked the radiation . the received energy values can be extrapolated to estimate the energy received at these locations outside the beam . one example is to set the received energy outside the right edge of the radiation beam equal to the received energy at the right edge of the radiation beam ; likewise , the received energy outside the left edge of the radiation beam can be set equal to the received energy at the left edge of the radiation beam . other extrapolation schemes are possible , however . finally , the reconstructed attenuation coefficient may be subject to a smoothing and de - blurring operation . the smoothing and de - blurring operation can be global , taking account of detected radiation across the whole field of view , or local , taking account only of detected radiation within a range of the location being considered . in the latter case , the requirement to extrapolate the received energy is reduced and therefore the accuracy of the reconstructed attenuation coefficient may be greater . all of these features are embodied in the equation ( 14 ) described below . however , embodiments of the present invention may implement only one or more of the algorithm features described above in order to reconstruct attenuation coefficients from the portal images . the present invention thus provides a method and an apparatus for reconstructing three - dimensional ct images from portal images acquired during treatment . such a possibility was previously thought impossible due to the narrow and ever - changing field of view , and the incomplete angular coverage of the reconstruction points by the rotating radiation source ( i . e . because of insufficient data ). however , embodiments of the present invention can be greatly simplified compared to the imaging techniques commonly employed conventional radiotherapy systems , which typically require a separate kilovoltage source of radiation and detector . there now follows a mathematical description of the algorithm employed in step 108 . the ‘ global ’ algorithm is shown in equation ( 10 ); the general algorithm , including the global and ‘ local ’ algorithms , is shown in equation ( 14 ). a hybrid version , designed to incorporate the benefits of both global and local algorithms , is shown in equation ( 16 ). there is denoted a 3d function describing the attenuation coefficient of a subject by the function , ƒ ( r , φ , z ), where ( r , φ , z ) are cylindrical polar coordinates . in practice , either the ram - lak kernel ( see ramachandran g n and lakshminarayanan a v 1971 “ three dimensional reconstruction from radiographs and electron micrographs : applications of convolutions instead of fourier transforms ” proceedings of the national academy of sciences us vol 68 , pp 2236 - 2240 , the contents of which are incorporated herein by reference ) or an apodizing kernel are applied in ct reconstruction . this results in the reconstruction not of ƒ ( r , φ , z ) itself , but rather this function smoothed by a point - spread function ( psf ). in the case of a ram - lak kernel this psf arises only from the finite pixel - width . in the case of an apodizing kernel , additional smoothing may be included to reduce image noise . consider such a reconstructed “ image ” in the ( r , φ ) plane . the smoothed 3d function , ƒ ( r , φ , z ), can be expressed as , where e r ( r ) is a 2d psf in the ( r , φ ) plane and ** represents a 2d convolution operation . now consider the cone - beam geometry illustrated in fig3 a and 3 b . the variables sad and sdd are the source - to - axis distance and source - to - detector distance , respectively . the variable β is the angle of the source with respect to the positive y - axis . the variables u and v denote the position of a point in the imaging plane . the feldkamp cone - beam reconstruction algorithm ( see feldkamp l a , davis l c and kress j w 1984 “ practical cone - beam algorithm ” journal of the optical society of america a vol 1 , pp 612 - 19 , the contents of which are incorporated herein by reference ) can then be written as , f r ⁡ ( r , ϕ , z ) = 1 2 ⁢ ∫ 0 2 ⁢ π ⁢ ⁢ ⅆ β ⁡ ( sad sad - s ) 2 ⁢ ∫ - u max u max ⁢ ⁢ ⅆ u ⁢ ⁢ d β ⁡ ( u , v ) ⁢ e r ⁡ ( u ′ - u ) ⁢ ( sad sad 2 + u 2 + v 2 ) , ( 2 ) u = u ⁡ ( sad sdd ) , v = v ⁡ ( sad sdd ) , ( 3 ) u ′ = t ⁢ sad sad - s ( 4 ) t = r ⁢ ⁢ cos ⁡ ( ϕ - β ) , s = r ⁢ ⁢ sin ⁡ ( ϕ - β ) ( 5 ) and d β ( u , v ) is the cone - beam ray - projection at the point on the detector defined by ( u , v ) for the source rotation , β . the function , e r ( u ), is related to the psf and takes the form , e r ⁡ ( u ) ≈ 1 2 ⁢ π ⁢ λɛ r ⁡ ( u ) ( 6 ) ɛ r ⁡ ( u ) = 1 r ⁢ π ⁢ γ ⁡ ( 5 2 + α ) γ ⁡ ( 2 + α ) ⁢ ( 1 - ( u r ) 2 ) α + 1 , if ⁢ | u | ≦ r = 0 , otherwise ( 7 ) where r defines a range and λ is the calderon operator . the calderon operator is pseudo - differential operator that performs a de - blurring operation . consider a euclidean space ″ in which a point is defined by the vector , r . a function in this space , ε ( u ), may have a fourier transform , e ( ρ ), and a fourier transform pair , e ⁡ ( ρ _ ) = f ⁢ { ɛ ⁡ ( u _ ) } = ∫ r ″ ⁢ ⅇ - 2 ⁢ π ⁢ ⁢ i ⁢ u _ · ρ _ ⁢ ⁢ ɛ ⁡ ( u _ ) ⁢ ⁢ and ⁢ ⁢ ɛ ⁡ ( u _ ) = f - 1 ⁢ { e ⁡ ( ρ _ ) } = ∫ r ″ ⁢ ⅇ + 2 ⁢ π ⁢ ⁢ i ⁢ u _ · ρ _ ⁢ e ⁡ ( ρ _ ) . ( 8 ) the calderon operator , acting on a function ε ( u ), is then defined as where λ acts on the n - dimensional space . in expression ( 6 ) the calderon operator acts on the scalar variable u only ( i . e . n = 1 ). note that it would also be possible to define e r ( u ) as a standard global ct convolution kernel , such as the ram - lak or hamming kernel . this would , however , preclude the extension to local ( lambda ) reconstruction discussed subsequently . returning to equation ( 2 ), it is noted that the factors appearing in the round - brackets of this expression are geometric weighting factors . it is assumed for simplicity that sad & gt ;& gt ; u , v and ignore such factors : this simplifies the resulting formulae . equation ( 2 ), however , is inadequate , as it stands , in dealing with the data insufficiencies in the projection data . the following vmat - ct formula is therefore proposed , f r ⁡ ( r , ϕ , z ) ≈ 1 2 ⁢ θ ⁡ ( ∫ 0 2 ⁢ π ⁢ m β ⁡ ( u ′ , v ) ⁢ ⁢ ⅆ β - β ′ ) 1 2 ⁢ π ⁢ ∫ 0 2 ⁢ π ⁢ m β ⁡ ( u ′ , v ) ⁢ ⁢ ⅆ β ⁢ ∫ 0 2 ⁢ π ⁢ [ m β ⁡ ( u ′ , v ) ⁢ ∫ - v max u max ⁢ ⁢ ⅆ ud β extrap ⁡ ( u , v ) ⁢ e r ⁡ ( u ′ - u ) ] ⁢ ⁢ ⅆ β , ( 10 ) where the extra constituents are defined in detail below . the mlc aperture is described by the masking function , m β ( u , v ). if the aperture extent , at a particular v , and β , is defined by the interval , u 1 ( v , β )& lt ; u & lt ; u 2 ( v , β ), then a possible form for the masking function is the acquired ray - projections should be extrapolated for the purposes of the convolution integral in ( 10 ). one possible extrapolation scheme is : d β extrap ⁡ ( u , v ) = d β ⁡ ( u , v ) , if ⁢ ⁢ u 1 ⁡ ( v , β ) ≦ u ≦ u 2 ⁡ ( v , β ) d β ⁡ ( u 1 , v ) , if ⁢ ⁢ u & lt ; u 1 ⁡ ( v , β ) d β ⁡ ( u 2 , v ) if ⁢ ⁢ u & gt ; u 2 ⁡ ( v , β ) . ( 12 ) more sophisticated extrapolation schemes are , of course , possible . this makes no assumption about specific mlc orientation . however , with some collimator rotations ( e . g . 90 °) it would be possible for leaves to split the aperture into more than a single region in the u - direction . such cases ( not considered here ) would require a slightly more complex masking function and extrapolation . the masking function is included inside the back - projection integral to prevent the back - projection of convolved ray - projections outside of the bev . however , to obtain imaging data of the target region , it is not important if the data set is not strictly complete over an angular range . that said , where the data is woefully insufficient the image quality will suffer , and thus a cut - off radian coverage , β ′, may be included . this parameter appears in the heaviside step - function , θ ( k ), of equation ( 10 ). the heaviside function forces the reconstructed value to zero at positions in the bev for less than β ′ radians . this allows the elimination of reconstruction points for which so little data are available that they are likely to be misleading . the denominator of the fraction in which the heaviside function appears is present to normalise the contributions according to the angular extent of data available . thus , for example , a voxel in the bev for 75 % of the 2π arc ( i . e . ( i . e . ⁢ ⁢ 1 2 ⁢ ⁢ π ⁢ ∫ 0 2 ⁢ ⁢ π ⁢ m β ⁡ ( u ′ , v ) ⁢ ⅆ β = 0 . 75 ) , ∫ 0 2π m β ( u ′, v ) dβ = 0 . 75 ), will be divided by a factor 0 . 75 to approximately correct for the missing contributions . due to data - insufficiencies in the set of projections obtained in a vmat treatment , exact and unique reconstruction of f r ( r , φ , z ) may not be possible . the data - insufficiency in the vmat - ct problem , in particular , the width truncation of projections by the mlc , suggests the possibility of using local tomography techniques ( see faridani a , ritman e l and smith k t 1992 “ local tomography ” siam journal of applied mathematics vol 52 , pp 459 - 484 , the contents of which are incorporated herein by reference ). the idea of local ( lambda ) ct is not to reconstruct ƒ r ( r , φ , z ) but rather a related object , where m is an integer denoting the number of repeat operations of the calderon operator . it is noted that in ( 13 ), λ acts on the n = 2 space of ( r , φ ) ( but not over z ). the new object ( 13 ) has many of the same properties as the original attenuation coefficient function . the generalised version of the vmat - ct reconstruction formula is then , λ m ⁢ f r ⁡ ( r , ϕ , z ) ≈ 1 2 ⁢ θ ⁡ ( ∫ 0 2 ⁢ π ⁢ m β ⁡ ( u ′ , v ) ⁢ ⁢ ⅆ β - β ′ ) 1 2 ⁢ π ⁢ ∫ 0 2 ⁢ π ⁢ m β ⁡ ( u ′ , v ) ⁢ ⁢ ⅆ β ⁢ ∫ 0 2 ⁢ π ⁢ [ m β ⁡ ( u ′ , v ) ⁢ ⁢ ⅆ u ⁢ ⁢ d β extrap ⁡ ( u , v ) ⁢ e r ( m ) ⁡ ( u ′ - u ) ] ⁢ ⅆ β , ( 14 ) e r ( m ) ⁡ ( u ) ≈ 1 2 ⁢ π ⁢ λ m + 1 ⁢ ɛ r ⁡ ( u ) . ( 15 ) if m is even , such as the m = 0 case , e r ( m ) ( u ) has a non - compact support and is termed “ global ”. the convolved ray - projection , to be back - projected , then depends on every unconvolved ray - projection along the u integration direction , for fixed values of v . therefore , if the patient does not entirely fit within the field - of - view at any angle β , required data are absent . despite this , in some cases , the ray - projections through the patient outside the acquired set can be extrapolated reliably enough to provide an acceptable reconstruction . it can be shown , however , that if m ≧− 1 and odd then e r ( m ) ( u ) can have a compact support and be described as “ local ”. this means that a convolved ray - projection , to be back - projected , depends only on the unconvolved ray - projections in its local vicinity ( within a distance r in fact ). thus the m = 1 algorithm has less restrictive data requirements than the m = 0 case . in some cases therefore we might expect local tomography ( m = 1 ) to provide more useful information than the global ( m = 0 ) algorithm . to encapsulate the benefits of either approach ( local or global tomography ) a hybrid semi - local construction is proposed , σ σ ƒ r ( r , φ , z )= ƒ r ( r , φ , z )+ σ rλƒ r ( r , φ , z ), ( 16 ) where σ is a dimensionless parameter and r is again the range parameter . it will of course be understood that many variations may be made to the above - described embodiment without departing from the scope of the present invention . although the present invention has been described with reference to preferred embodiments , workers skilled in the art will recognize that changes may be made in form and detail without departing from the spirit and scope of the invention .
a radiotherapy system comprising a support for a patient undergoing radiotherapy treatment , a gantry rotatable about an axis , a source of radiation mounted on the gantry and producing a beam of radiation directed towards a target region of the patient , a collimator coupled to said radiation source , the collimator comprising a plurality of movable , beam - limiting elements , to collectively define a shaped aperture through which the radiation beam passes , a portal imager mounted on the gantry opposite the radiation source for detecting the radiation after it has passed through the patient and generating corresponding images , and associated circuitry for controlling at least the gantry , the source , the collimator , and the portal imager , collating detected data comprising a plurality of images acquired including images at a plurality of angles of rotation of said gantry and images at a plurality of collimator shapes ; generating a three - dimensional image of the target region based thereon .