Method of Reconstructing a High Resolution Image by Tomodensitometry

1. Field of the Invention

The invention relates to a method of reconstructing a high resolution image by tomodensitometry with the purpose of improving the quality of the reconstructed image without burdening the technological structure of the means for measuring the penetrating radiation absorption and without substantially modifying the data acquisition procedure.

2. Description of the Prior Art

So called "third generation" tomodensitometers are essentially formed from a source of penetrating radiation (X rays) emitting a flat fan-shaped beam, a curved row of detectors placed opposite this source so as to measure the unabsorbed fraction of the radiation in the plane of the section of which an image is to be formed, means for rotating the source and the detectors in the plane of this section and a processor for processing the data supplied by the detectors and for reconstructing an image from this data. The spatial resolution of the reconstructed image depends more particularly on the number of detectors covering the fan shape of the X ray beam. For a given beam aperture, the higher the number of detectors the better the spatial resolution. However, the multi-detector assembly is one of the most expensive parts of the tomodensitometer and one of those which is the most difficult to produce with good reliability. Consequently, the increase in the number of detectors on a given sector corresponding to the fan shape, comes up against technological limits.

One of the aims of the invention is to substantially improve the quality of the image without modifying the technological structure of the source-detectors assembly and without substantially modifying the data acquisition procedure properly speaking.

Another aim of the invention is also to reduce the cost of construction of a tomodensitometer having spatial resolution comparable to that of existing tomodensitometers, this reduction in cost being obtained by using a multidetector assembly comprising only half the number of detectors of exising machines.

For this, the invention provides essentially a method for reconstructing images by tomodensitometry, of the type consisting in rotating a source-detectors assembly in the plane of a section the image of which it is desired to form, said source emitting a penetration radiation beam in the form of a fan and taking and storing in a memory a succession of views for different rotational angles .beta. of the source-detectors assembly with respect to a fixed reference axis system of said plane, the detectors being adapted so as to converge towards the focus of said source and each view being sampled for angles .gamma. representative of the detectors and referenced from the straight line joining the center of rotation of the source-detectors assembly and said focus, the angular shift between two adjacent detectors .DELTA..gamma. being constant, and the angular shift between views .DELTA..beta. being constant, which method consists, with the angular shift between said straight line and the nearest detector thereto being adjusted to ##EQU2##

in taking and storing a series of views spaced apart over a complete revolution of said source-detector assembly;

in working out for each view intermediate values representing an interlaced sampling with that of said detectors, each intermediate value, of a given angle .gamma., being obtained by a combination of at least two linear absorption values inferred from measurement values taken from views situated on each side of a fictitious view whose image by transformation

contains the point representative of said intermediate value in space [.theta.,.gamma.] with .theta.=.beta.+.gamma. and k a relative integer, said values of measurement of said views being taken from the angle detector opposite the one of said desired intermediate value; and

in applying a reconstruction algorithm, known per se, comprising a convolution and back projection to the views thus completed by said intermediate values.

The invention will be better understood from the following description of a method based on its principle and of the essential component parts of a tomodensitometer operating in accordance with this method, with reference to the accompanying drawings in which:

FIG. 1 is a schematical representation of the main component parts of a tomodensitometer for implementing the invention;

FIG. 2 is a representation in space [.theta.,.gamma.] of the views taken by a conventional tomodensitometer operating in the pulsed mode for reconstructing an image;

FIG. 3 is a representation in the same space [.theta.,.gamma.] of the views taken by a tomodensitometer of the same type in which one of the modifications for implementing the invention has been made, namely a static shift between the source and the detectors;

FIG. 4 illustrates, in the space [.theta.,.gamma.]

FIG. 4 illustrates, in the space [.theta.,.gamma.] a transformation which may be applied to the views of FIG. 3 for implementing the method of the invention;

FIG. 5 illustrates in the space [.theta.,.gamma.] the result of this transformation and shows the interpolation procedure for completing the views;

FIG. 6 illustrates the use of an algorithm in accordance with the invention, for calculating an intermediate value of the first view V.sub.1 ;

FIG. 7 illustrates a phase in the calculation of the same algorithm; and

FIG. 8 is a representation in the space [.theta.,.gamma.] of a view in the case of a tomodensitometer operating in continuous mode.

In the rest of the description with reference to FIGS. 1 to 7 the particular case is taken of a tomodensitometer operating in pulse mode. This kind of tomodensitometer is well known. In most cases, the pulsed mode is given concrete form by the fact that the fan shaped beam is only emitted for a very brief period of time at the moment when the source-detectors assembly is at a predetermined angle .beta.. In this case, it is sufficient for the reading of all the detectors to be finished when the source-detectors assembly reaches the next angle .beta.. However, another way of providing said pulsed mode may consist in permanently generating the fan shaped beam but in "reading" all the detectors simultaneously at each angle .beta..

Referring to FIG. 1, the essential component parts of a tomodensitometer have been shown. This apparatus comprises a source-detectors assembly rotatable in a sectional plane, the image of which it is desired to form. The center of rotation is 0 and the plane of the section will be considered as being that of FIG. 1. With this sectional plane is associated a fixed reference orthonormed axis system XOY. The source is shown symbolically in the drawing by its focal plane F. It is, as is well known, a source of penetrating radiation and more particularly an x ray source emitting a fan shaped beam 11. In the rest of the text, mention will be simply made of source F. The case 13 containing all the detectors d.sub.i has a curved shaped so that these detectors are arranged along an arc of a circle CC' centered geometrically on the focal point F. The case may of course be rectilinear but will then contain detectors arranged so as to converge towards the focal point. In the example described, the fan shaped beam 11 covers an arc of a circle CC' occupied by the input faces of the detectors d.sub.i. Thus, the reconstituted image will be inside circle 16 with center 0, inscribed within the limits of angle CFC'. An assembly may be used comprising 1024 detectors d.sub.i. Case 13 and source F are fixed to the same rotary support (center of rotation 0) and are stationary with respect to each other so that rotation of the source-detectors assembly may be referenced by the angle .beta. between axis OY and a straight line OF. The attenuation values of the x ray beam for predetermined angles .beta. are inferred from the values of measurements taken at the outputs of the detectors and stored in a memory M of a processor. For each predetermined angle .beta. an order for reading the detectors d.sub.i is emitted, for example by means of an incremental coder 16 adapted for reading angle information related to the rotating source-detectors assembly and for delivering (drive connection 17) orders for controlling the reading of the status of all the detectors. The whole of the linear attenuation values inferred from reading the detectors for a given angle .beta. is called "view".

Thus, it may be considered that each view is sampled by the very structure of the detector assembly, for angles .gamma. representative of the detectors. Each angle .gamma. indicating the position of a detector d.sub.i is referenced from the straight line OF. Thus, the constant angular shift between the centers of the input faces of two adjacent detectors will be called .DELTA..gamma..

If, for each angle .beta. of a view, we consider the axis OX perpendicular at P to the ray of angle .gamma. striking any detector, it can easily be shown that the angle .theta. between the axis OX of the fixed reference and the axis OX is such that:

and this for all cases, since the angles .theta., .beta. and .gamma. are orientated. We have again:

From this latter relationship, it can be seen that the different views (for different angles .beta.) have an extremely simple representation in space [.theta.,.gamma.].

As shown in FIG. 2, these views are carried by straight lines parallel to the first bisectrix and shifted by .DELTA..beta. along the axis of the .theta.. Points may be plotted on these straight lines representing the coordinates in the space [.theta.,.gamma.] of the sampled values effectively acquired on reading the different views.

Thus, assuming that the number of detectors is even, i.e. 2n and that the straight line OF passes exactly between the two center detectors d.sub.-1 and d.sub.0, they appear shifted by ##EQU3## on each side of the axis .gamma.. The aim of the invention is to obtain for each view intermediate sampled values whose coordinates in the space [.theta.,.gamma.] are represented by stars * in FIG. 2. In other words the aim of the invention is to obtain views comprising 4n sampled values with an assembly of 2n detectors only and thus to improve the spatial resolution of the reconstructed image.

According to a feature of the invention, a permanent and predetermined angular shift between the assembly of detectors (i.e. case 13) and straight line OF is caused above all. In the example chosen, the angular shift between the straight line OF and the nearest detector (d.sub.o) is chosen equal to +(.DELTA..gamma./4). That means that the assembly of detectors is shifted "to the left" considering FIG. 1, so that the center of the input face of detector d.sub.o is shifted angularly by +(.DELTA..gamma./4) with respect to OF and so that the detector d.sub.-1 is shifted by -(3.DELTA..gamma./4) with respect to the same straight line.

FIG. 3 shows the configuration of the views in the space [.theta.,.gamma.] after this shift.

For implementing the invention, a series V.sub.i is taken of N views spaced apart over a complete revolution of the source-detectors assembly, with N=(2.pi./.DELTA..beta.).

The desired intermediate values s.sub.i are also located on the different views in FIG. 3, each intermediate value position in the space [.theta.,.gamma.] being equidistant from the positions of two adjacent real detectors.

That is to say the transformation:

It can be shown that the function of absorption of a view is invariant for this transformation T. Now, this transformation T may be broken down in the space [.theta.,.gamma.] in to translation -.pi..theta. followed by a symmetry with respect to .theta.. This transformation is illustrated as such in FIG. 4. If said transformation is applied to the view V.sub.i, the result is the "view" V'.sub.i.

If the transformation T is thus applied to the series off views V.sub.i, it can be seen that the sampled values of the transformed "views" V'.sub.i have the same ordinate .gamma. as the desired intermediate values s.sub.i. Thus the views V.sub.i may be completed by interpolations from linear absorption values corresponding to sampled values of the "views" V'.sub.i, with .gamma. constant, in the space [.theta.,.gamma.], as shown by the arrows in FIG. 5. In the example, each intermediate value of a given angle .gamma. is obtained by a combination of two linear absorption values inferred from the sampled values having the same ordinate .gamma., taken from the "views" V'.sub.i situated on each side of the desired intermediate value. The combination may simply consist of a linear interpolation between the corresponding linear absorption values of two adjacent "views" V'.sub.i, as shown in FIG. 5. In other words, each intermediate value is obtained by averaging these two above mentioned linear absorption values (s.sub.o ' and s.sub.o "), by weighting them respectively by the inverse (1/W' and 1/W") of the difference (in projection along the axis of the .theta. FIG. 5) between their respective positions (s.sub.o ' and s.sub.o ") and the position of the desired intermediate value (S.sub.o). For some "views" V'.sub.i required for the interpolation calculations suggested by FIG. 5, transformation T uses a view V.sub.i situated outside the range [.theta.,2.pi.] of the views actually acquired. However, these latter may be theoretically extended in the space [.theta.,.gamma.] from -.infin. to +.infin. because of the periodicity, with respect to .theta., of the linear attenuation function beyond a complete revolution (because of the condition N=(2.pi./.DELTA..beta.), which may be represented thus:

In other words, the transformation T may be completed in the following way:

with k a relative integer.

There always exists a value of k such that the view V.sub.i on which the translation .theta.+(2k-1).pi. is applied is indeed one of the views acquired V.sub.i. FIG. 5 illustrates the justification of the procedure of the invention. This method of calculating the intermediate values must now be explained for the "views" V'.sub.i only represent a rearrangement of the views V.sub.i and the calculation of each intermediate value may in fact be made from each "position" of the desired intermediate value in the space [.theta.,.gamma.]. This is what is shown in FIG. 6 where only the first view V.sub.1 and two "views" V'.sub.i have been shown, a result of the transformation T applied to the views V(N/2) and V(N/2)+1 (N being even), which are also represented. The views V(N/2) and V(N/2)+1 are situated on each side of a fictitious view VF.sub.o whose transform VF'.sub.o, by the transformation T, passes through the coordinates of the intermediate value s.sub.o of the view V.sub.1 in the space [.theta.,.gamma.]. According to the invention, it is therefore sufficient, for each intermediate value (S.sub.o) of each view (V.sub.1), to search for the corresponding "fictitious view" (V.sub.Fo) whose transform by the transformation T passes through the coordinates of said intermediate value, to select the two adjacent real views V.sub.i, V.sub.i+1 (V.sub.N/2, V.sub.N/2+1) situated on each side of said fictitious view and to take from these two views the absorption values of angle - .gamma.(.gamma. being the ordinate of s.sub.o in the space [.theta.,.gamma.]) which correspond to measurement values from the same real detector and to effect a linear interpolation between these linear absorption values. It should be noted that, for the sake of simplicity, it was assumed above that the number N of acquired views V.sub.i was even; the procedure of the invention can nevertheless be applied such as it is if N is uneven.

One of the possible calculation algorithms will now be described, forming a presently preferred embodiment of the above described method. In the example described, we have the case of a "leftward" shift of the detectors, such as illustrated in the drawings and the case of a pulsed operating mode where, as is known, the sampled values of the views V.sub.i are carried by straight lines parallel to the first bisectrix of the space [.theta.,.gamma.]. With reference to FIG. 7 the algorithm will be first of all described for completing the view V.sub.1. The equation of view V.sub.1 is:

It is desired to calculate the intermediate value as s.sub.j.

Let x.sub.j be the difference between the intersection of the view to be completed with axis .theta. (in this case the origin, in the case of view V.sub.1) and the interconnection of the transform of the fictitious view (passing through s.sub.j) with the same axis .theta..

We have:

The point of intersection of the corresponding fictitious view with the axis of the .theta. is therefore:

The number of the real view KG situated "on the left" of this fictitious view is therefore determined by the whole value of the ratio (y.sub.j /.DELTA..beta.).

More precisely: ##EQU4## and the number of the view KD situated "on the right" of the fictitious view is therefore:

It now remains to determine "the ordinate" .gamma.d.sub.i of a detector .gamma..sub.i to be taken into account in each view KG and KD, i.e. according to the transformation T:

Now, if we refer to FIG. 3:

and

i.e.

Further, let r.sub.j be the rest of the operation ##EQU5## r.sub.j represents the difference in projection along the axis of the .theta. between the position of s.sub.j and that of detector d.sub.j in the real view KG; r.sub.j is shown in FIG. 7 between the transform of the fictitious view passing through s.sub.j and the transformed KG' of KG. Consequently, from r.sub.j can be inferred the weighting to give to the linear attenuation values corresponding to the measurement values taken from the detector d.sub.i of the views KG and KD.

Let h.sub.KD (d.sub.i) and h.sub.KG (d.sub.i) be the linear attenuation values corresponding to detector d.sub.i in the real views KD.sub.(j) and KG.sub.(j), and let:

It can be seen that the linear attenuation value h(s.sub.j) to be attributed to the position s.sub.j of view V.sub.1 (in other words the intermediate value s.sub.j itself) is expressed:

Programming the above expressed A to F is within the scope of a man skilled in the art, as is also the design of a specialized calculating processor coupled to the buffer memory M.

To complete the following views V.sub.m (i.e. V.sub.2, V.sub.3, . . . , V.sub.N) we have, instead of the relationship (A):

that is to say:

The expression A' replaces A in the rearrangement and interpolation algorithm but the other relationships B to F remain unchanged. It is also interesting to note that from one view to another R.sub.(j) remains unchanged and that KG.sub.(j) is incremented by 1. Consequently, it is possible in the computing software to precompute r and KG solely for the first view and then to procede by incrementing KG. Furthermore, it can be seen that the same algorithm is valid for a "rightward" shift of the detectors, so that the method described can be applied for a shift of the detector assembly of: ##EQU6##

The preceding algorithm was developed on the assumption of pulse mode operation. It is easy to adapt it to continuous mode operation. According to a usual continuous operating mode, it is assumed that source S emits permanently and that reading of the measurement values at the outputs of the detectors, for each view, is effected by progressive scannng from one end to the other of the assembly of detectors, for example from detector d.sub.-n to detector d.sub.n-1.

Two additional parameters must be taken into consideration:

the speed of rotation v of the source-detectors assembly;

the time interval .tau. separating the acquisition of the sampled values corresponding to two adjacent detectors, assuming a constant scanning speed.

FIG. 8 illustrates the way in which the representation of a view in the space [.theta.,.gamma.] evolves if the parameters v and .tau. are taken into account. The view carried by the broken straight line corresponds to the pulsed mode whereas the same view carried by the continuous straight line corresponds to the continuous mode.

The sampled values move parallel to the axis .theta. since .theta.=.beta.+.gamma.; .beta. varying during the acquisition. The two straight lines converge at d.sub.-n since scanning of the detectors begins by this detector and the next sampled value d.sub.-n+1, is shifted at d'.sub.-n+1 by a distance v..tau. in the space [.theta.,.gamma.], parallel to the axis of the .theta.. Let a be the angle between d.sub.-n d.sub.-n+1 and d.sub.-n d'.sub.-n+1. ##EQU7##

The views are therefore carried by straight lines parallel to the straight line of equation .theta.=(1+.alpha.).gamma.

Consequently, for the view V.sub.1, the relationship (A) becomes:

and for any view V.sub.m, the relationship (A') becomes:

The rest of the algorithm is applied without modification with the new above defined values of x.sub.j.

It can be seen that the method described can also be applied without change whatever the type of continuous mode scanning. In this case, the sample values of the views are no longer necessarily carried by straight lines but by any curves to which the same transformation T may be applied so as to infer therefrom the desired intermediate values from interpolations similar to those described above.

Referring again to FIG. 1, the essential parts of the computer associated with detectors d.sub.i are shown in the form of a functional block diagram. The sampled values read at the outputs of the detectors d.sub.i are transformed into digital information by means of an analog-digital converter 18 then processed in a preprocessor 19 carrying out the usual calibration and logarithmic transformation operations. The information thus transformed is formed by the corresponding linear attenuation values and this digital information is stored in memory M.

For reasons which will be more clearly understood further on, memory M is in the form of a matrix of storage units comprising 4n columns and N lines (number of views V.sub.i) since the acquisition of the views V.sub.i is made during a complete revolution and with: .DELTA..beta.=(2.pi./N).

Filling of memory M with the information representative of the views V.sub.i occurs on the basis of one storage unit out of two; these units are shown by hatching in FIG. 1. Filling of memory M, such as shown, corresponds to a "leftward" shift of the detectors. Memory M is associated with reading means 20 which organize the transfer of information to a processor 21, known per se, and whose role is to apply the above mentioned reconstruction algorithm. However, reading of a line of memory M is only ordered when all the storage units of this line which have not been "filled" during the phase for storing the views V.sub.i, have received digital information representative of the intermediate values, calculated using the above described algorithm. For that, a specialized processor or a sub program of the computer effects the computing, addressing in memory M and linear interpolation operations before ordering writing of the resulting intermediate values in the storage units left temporarily "empty" during writing of the real views in the same memory M.

This processor or this sub program is shown symbolically in FIG. 1 by the functional block 22 adapted for calculating the addresses of the useful linear attenuation values, reading these values, effecting the interpolation calculation combining these values and rewriting the result at the suitable address of an intermediate value of a view to be completed. It also drives the reading means 20 (functional connection 23) controlling the transfer of information from memory M to the reconstruction processor 21, the transfer of a "line" from the memory only being permitted when all the storage units corresponding to the intermediate values have been filled.