Method of correcting the geometric distortion of an image, in particular of an image produced by a gamma camera.

The present invention relates to a geometric distortion correction method of an image, especially an image produced by a gamma camera. It finds its application particularly in the medical field where gamma cameras are used as a diagnostic aid. It relates to scintillation cameras (or gamma camera) type ANGER including US patent 3,011,057 describes the operation in its principles and means of implementation. These gamma cameras are designed to detect and visualize the photons emitted by the radioactive body.

Gamma cameras are used in nuclear medicine to visualize an organ in the distribution of molecules labeled with a radioactive isotope which has been injected into a patient. A gamma camera generally comprises a collimator to focus the gamma photons emitted by the patient, a scintillating crystal to convert the gamma photons into light photons or scintillations and a network of photomultiplier tubes which convert each scintillation into electrical pulses said electrical contributions tubes. They also contain electronic circuitry to produce from electrical inputs from the photomultiplier tubes, X and Y coordinate signals from the place where occurred the scintillation and a validation signal Z when the energy W of the scintillation belongs to a predetermined energy band.

This detection channel is generally followed by a display assembly generally comprising a cathode-ray oscilloscope controlled by the X signal, Y, and Z for viewing by a luminous point on the point of impact of the gamma photon on screen crystal. This impact is also called image event. The display unit may optionally comprise a photographic apparatus for forming an image of the observed body by integrating a large number of light spots produced on the CRT screen. It may also include a digital processing device images. In particular the display unit can be adapted to the presentation of tomograms of the organ observed. To achieve this goal we acquire multiple images of the organ at a plurality of observation of the gamma camera orientations relative to this body. By signal processing, similar to those encountered in CT, we can reconstruct section images of the examined organs.

Among other qualities a gamma camera must have good resolution Spatials, ie the ability to distinguish small radioactive sources close, a good response in counting rates, ie the ability to process a large number of of events per unit time, and an independent image quality of the energy of the isotope in question. The spatial resolution depends on the accuracy of the calculation of the X and Y coordinates of each of the image events. The quality of the development of these coordinates depends mainly on the physical laws governing the operation of various parts of the gamma camera. Thus the interaction of a gamma photon with the crystal gives rise to a light scintillation whose intensity decreases exponentially with time. The time constant of this decay is characteristic of the scintillator crystal used. For example, for a sodium iodide crystal thallium-activated, NaI (T1), it is of the order of 250 nanoseconds. This scintillation is seen by several photomultiplier tubes simultaneously. The light photons forming this scintillation tear photoelectrons from the photocathodes of the photomultiplier tubes. The number of photoelectrons obeyed for a given scintillation, statistics Poisson. This means that the electric contribution of a photomultiplier tube receiving a scintillation has an amplitude whose value follows a statistical Poisson distribution and the value of which depends on the energy of the light incident photons. Moreover, at constant energy, the electrical contribution is substantially a Gaussian function of the distance separating the center of the photomultiplier tube of the place where is the scintillation occurs. If the scintillation takes place in line with the center of this tube, the electrical contribution is maximum. More venue of the scintillation is from the center of the tube, the more power contribution is low. For example, if a scintillation occurs in vertical alignment with a wall of tube the electrical contribution thereof is approximately halved with respect to the maximum electrical contribution.

A scintillation is seen by several photomultiplier tubes simultaneously, usually six to ten tubes. Also, determining the location of said scintillation on the crystal itself represents the source of emission of the gamma photon excitation (and thus the image event) may be obtained by calculating the location the barycenter of the electrical contributions supplied by all the photomultiplier tubes excited by said scintillation. This calculation is done simply by ANGER, by injecting the electrical contributions through a set of resistors of matrices whose resistance values ​​are a function of the positions of the photomultiplier tubes to which they are connected. The positions of these tubes are located with respect to reference Cartesian axes whose point of intersection is generally centrally located in the network of tubes. In each matrix, there are as many resistors as there are photomultiplier tubes in the tube array. Each of the resistors is connected firstly to the output of a different photomultiplier tube and the other to a common point which forms the output of the matrix. These resistors and carry a weighted electrical contributions of each of the photomultiplier tubes supplying them.

One of the problems presented by the detector gamma cameras of the Anger type is that they have geometric image distortions related to the light-gathering structure: crystal scintillator - photomultiplier tubes - barycentration matrices. Advances in nuclear medicine, particularly the improvements made ​​to gather more information and better information from gamma cameras, eg for the detection of small tumors, leading to a lack of linearity inherent in the design and the achievement of cameras. This results in a spatial distortion of disturbing images. These distortions can cause significant density uniformity defects. For example a contraction of 0.4 mm radius of the image of a surface of circ ular radius of 1cm leads to a 8% non-uniformity. These non-uniformities can be particularly troublesome in tomography, where the amplification effect of the nonuniformities due to the methods of reconstruction may be greater than a factor of 10.

We tried in the prior art, to correct the effects of this distortion. The general principle of the correction is as follows. Measurement is carried out of an image of a regular pattern disposed between a radioactive source of uniform emission and the gamma camera. If the gamma camera was perfect, follow the picture should match the regular distribution of holes or slots of the target. The measurement actually shows irregularities from the spatial distortion effects. However, one can use in the image thus acquired knowledge of the distortion, measured by comparing the acquired image to the theoretical distribution to which should have led, to perform corrections subsequently acquired images with the same gamma camera. Thus, in US Patent 3745345, the test pattern, the ghost of acquisition comprise a series of 3 mm diameter holes spaced each other of 24 mm. Four acquisitions were made by shifting each time X and Y, the target of 12 millimeters. For each measuring point a pair of correction coefficients is calculated and stored. Intermediate points are determined by interpolation up to 64 x 64 for the entire image field. In normal operation, the distorted image is acquired and stored. At the end of acquisition, a correction program redistributes the image Events distorted image array in a corrected image matrix obtained through correction factors. This method has two drawbacks. On the one hand the redistribution of the image events depending on the coated surface leads to enlargement of the content of a mesh on four meshes. This automatically results in a loss of spatial resolution of the corrected image: it becomes less accurate. On the other hand, because of the distortion, the number of blows to redistribute is not accurate. Indeed, there is no reason why the surface of a unit cell in the acquired image data and distorted is equivalent to the surface of a unit cell in the corrected image.

In another state of the art represented by the French patent 2,412,856, the ghost of acquisition is a 3 mm slots pattern spaced 15 millimeters. Four acquisitions were also made by shifting 7.5 mm in X and Y. Each X and Y coordinates of the distorted image, a program calculates the new coordinates U and V of the undistorted image. This program does not calculate a correction factor directly but the new coordinates. These new U and V values ​​are stored in a 64 x 64 matrix In operation correction, X and Y are coded on 12 bits. The six bits of X and Y memory address U and V. Then a linear interpolation is carried out using the six least significant bits of X and Y is thus used the shift of an image event to place over the four corners of the unit cell of the corrected image in which the image event must ultimately be. Given a picture field about 400 millimeters in diameter, a codification of twelve bits of image coordinates events should normally lead to a picture of precision in the order of 0.1 millimeters. However, experience shows that it does not lead to the accuracy of the corrected image. The document calls 2412856 well replace linear interpolation by non-linear interpolation, for example, a cubic type, but this document adds that though the linear interpolation provides sufficient accuracy to calculate the true coordinates. This means in fact, and the experiment verifies that the non-linear interpolation does not give better results. We shall see later, the expected accuracy of 0.1 millimeters in the erect image, thus can not be reached.

A third prior art, consisting of the European Patent 0021.366, recommends using a phantom acquisition, which is a pattern of slits of 1 mm, spaced 15 millimeters. Data is stored in a matrix of 256 x 256 with 12 bits. A computer produces a matrix of 64 x 64 correction correction coefficients. Again, six bits memory address correction coefficients and allow to succeed in elementary mesh correction which are also associated eight coefficients corresponding to a finer correction of interpolation which uses six bits low weight of the calculated coordinates. Subsequently, a consistency analysis is done in phase calibration and correction. The computer calculates the events Fi image density per unit area Ai of the distorted image acquired. Then it calculates a F'i density corresponding to an elementary surface of the corrected image of A'i surface obtained after application of the correction coefficients. If distortion correction was perfect F'i should be uniform. The cited European patent recommends using the gradient of the corrected density for iteratively modify the coefficients of distortion correction. This patent thus highlights what appeared to reading the previous two, that it is essential to consider the information density depending on the distortion. The solution recommended in the latter document also presents a significant drawback in that it introduces defect density of the image, inherent in other reasons as image distortion, in correcting image distortion. Once applied this technique, it follows that the expected accuracy can not be attained because the improvement tends to be spread over several adjacent meshes of the defects that arise as a mesh, without otherwise discriminate reasons.

the question of why these techniques are then arose which, given the number of bits to encode utlise the coordinates of the image events would normally result in an expected image detail are not resulted. In particular, we posed the question of why, in the second state of the cited art, the correction approximation by cubic polynomial was no better than linear interpolation. then it was discovered that the choice of mesh size 64 x 64 for a whole image field diameter of 400 mm is defined in a too large mesh which provided thereby have an interpolation basis from the start of insufficient sampling the distortion function and could deduct the number of bits which were coded coordinates of the image events. On the other hand it was appercu as interim interpolation calculations devainet be made with greater precision of 1/10 of a millimeter if you wanted at the end of calculating the desired precision, ie 1 / 10 ° of a millimeter to a residual non-uniformity of the order of 2%. Ultimately, it was determined that there was a relationship between the minimum number of unit cells in the mesh used to establish the image correction matrix and the expected accuracy of the spatial distortion correction; this relationship is a function of spatial distortion gradient. It became so that it was necessary to make a more accurate interpolation calculation the required minimum accuracy.

Also the invention relates to a spatial distortion correction method of the image produced by a gamma camera, this gamma camera providing output signals representing, in an image control, image event coordinates , the method comprising the steps of:     - Measuring the coordinates of the image events obtained by interposing a test pattern corresponding to a theoretical uniform distribution of these images events     - By cubic approximation is calculated the coordinate displacements of the points of a theoretical mesh which corresponds to the test image,     - Deduced distortion corrections to be assigned to the coordinates of the image events,     - And applied to these coordinates, characterized in that     - Determining a minimum number of unit cells in the mesh corresponding to the image field, according to a required accuracy of the correction of spatial distortion.

The invention will be better understood from reading the following description and examining the accompanying figures. These are indicative only and in no way limit the invention. The figures show:<ul><li>- 1: a test pattern used to implement the method of the invention;</li><li>- Figure 2: Profile of the detected signals for one line from the image of the pattern;</li><li>- Figure 3: comparing a parallel lines equidistant network corresponding to the theoretical image to be obtained at the distorted image actually obtained;</li><li>- Figure 4: calculating by cubic approximation of the coordinate displacements of the points of a mesh of the test image;</li><li>- Figure 5 is a particular method of calculation of the unit cell according to the invention;</li><li>- 6: the interpolation method used in the invention to correct an image</li></ul>

1 shows a test pattern to be used to implement the invention. The test object 1 has slots such as 2, whose width is of the order of a millimeter and which are spaced from each other by strips 3, about twelve millimeters wide. Given an image field diameter of 400 millimeters, there are thirty three slots. The center slot is slightly offset from the center of the target of 3 mm so that if she returned to 180 degrees, we can not define another interlaced with the first slits. It is assumed that this pattern is placed between a not shown gamma camera (equal image field) and a density gamma ray radiation source emissions nearly uniform in the image field. The image of this pattern is analyzed according to a number of lines. For example the ordinate Yi of a line, one observes, at the end of a certain duration, the histogram H of the events of images based on the X abscissa of the points of the line.

The histogram H shown in Figure 2 has the appearance of a succession of peaks approximately spaced further step 3 of the pattern 1. If m is called the maximum of these histogram peaks, we note that for frequencies 0.35 m the width of the ridges is about four millimeters. It becomes unnecessary to search result to the image of a pattern of finer pitch: runs the risk that interpenetrate ridges on adjacent slots. For each of these peaks is calculated the abscissa Xc of center of gravity of the peak. Given that was used thirty three slots thus obtain thirty three values ​​of Xc. After reversal of the pattern, we can have thirty three other Xc values ​​in principle thirty interlaced with the previous three. Calculating the abscissa of centroids is of a known type. Preferably it is carried out by digital processing program in a computer.

thus calculated for each of Yi ordered lines of the test image 66 abscissas of centroids. In practice, this calculation to 512 lines using the 9 bits of the ordinates of the image events considered. 3 shows in the form of series of points alignments of these centroids according to their Yi ordinate line. Obviously they correspond to slots, the central axes a1 a2 etc ... are equidistant from each other. The determination of the axis positions have relative to the centroid X c Yi is performed taking into account the equidistance one hand and applying a method of least squares on the other. According to this method the sum of squares of DX abscissa variations is minimized for a set of axes have parallel and equidistant from each other. The program that performs determining the network straight've done over a line by line analysis to see if there is no lack of centroid. If any are missing, they are determined by linear interpolation between the neighboring centers of gravity corresponding to the same slot have.

The program then performs a horizontal interpolation by a set of cubic polynomials in order to reach 512 points from the 66 abscissa centroids measured Xc. Figure 4 shows schematically this interpolation and shows the effect of the invention. This figure represents the displacement DX abscissa centroids to Xci place abscissa calculated centroids of these. By a set of cubic polynomials one is able to reconstruct the curve 4 on which are determined the values ​​DX 512 now 512 corresponding to values ​​of Xc. In contrast to the state of the art cited, it was discovered in the invention that to achieve the expected accuracy (in the case of twelve bits and 0.1 mm) was necessary, given the order magnitude of the distortion gradient generally obtained on camera ANGER, using a unit cell of the correction mesh that is of the order of 1/128 of the image field. In the prior art cited a mesh 64 X 64 was used. As we shall see later, a 64 x 64 mesh led to significant residual defects in image areas with high distortion gradient. This mesh and intuitively corresponds approximately to the number of the test pattern used slots. The cubic approximation the writing was done in the prior art tended to replace the abscissa measures barycenter, for example 66 of the example mentioned, by sixty four correction abscissa which for some were relatively close measured values ​​as 66 is very different from 64.

Then, to complete the entire correction in the prior art cited, was used in a unit cell thus defined, linear interpolation. But this meant that the corrected abscissa, obtained by linear interpolation, were tainted from the outset of an error that in some cases was zero (for example at the location of point 5) and in other cases was very important (e.g. error 6 at the location of the abscissa 7). This error may be greater than the expected accuracy, the interpolation calculation performed later was illusory because it was on the wrong foot. Against by, in the invention, one chooses a finer unit cell, 128 X 128, so that the curve 4 is precisely approximated, and in such a way that linear interpolations on a finer mesh lead to structural errors 8 below the expected precision. In this case it was found that it was necessary to choose a mesh size equal to or finer than 128 X 128.

To perform the cubic approximation equations the limits are such that the cubic polynomial passes through 9, 10 and 11, and 12 that the derivative of the polynomial is the bisector 10 of the segments 9 - 10 and 10 - 11. The method of calculating of abscissa corrections by a cubic polynomial provides DX values ​​as a function of Yi, DX (Yi), which follow fairly continuously in function of the abscissa X. against, Yi of a line to another line immediately adjacent the DX values ​​to the same abscissa are relatively discontinuous. The program can then perform a vertical smoothing coefficient variable according to the variation of the slope in the vertical direction. On the edges of the image, the slots are extrapolated by parables. Moreover, the very shortened slits at the edges of the field, right or left, are extrapolated by straight. This results in a marked improvement in distortion correction on the edges of the image. This is not very useful in direct imaging (the practitioners always manage to put the body to visualize the middle of the image field). However, this is very useful in tomography extended organs where it actually reduces the artifacts rings.

At that stage, only half the work was done. It is then necessary to turn the pattern 1 of 90 ° (then again to 180 °) and repeat the procedure for calculating the corrections ordered DY (Xi) according to Y. This work is done as above for X is finally obtained two correction coefficient matrices of 512 X 512. given the foregoing it returns to a mesh size of 128 X 128 by averaging, as shown in Figure 5, the abscissa corrections to be applied to each of the sixteen elements of adjacent images of the mesh 512 and X 512 belonging to the same mesh 13 mesh 128 X 128.

When an image is acquired can decide to do with the correction matrix correction coefficients 512 without additional interpolation. The image uniformity obtained has a moire which is due to recovery or to the space between the neighboring pixels as there is no interpolation. This it is necessary to perform an interpolation to avoid moiré. This interpolation principle for the identification of a unit cell of the distorted image in which there is an image event to be seen in the corrected image. Let ABCD, 6, the four points in the corners of the distorted mesh elementary question. The correction coefficients for defining the corrected elementary mesh materialized in Figure 6 by the dashed segments 14 to 17. If we call DX DY and the two correction coefficients applied to each of the points A to D and concerning respectively segments 14 to 17 can be calculated, knowing the xy coordinates of the image P event in the unit cell ABCD image distortion corrections, DX and DY, to be applied to XY coordinates of the point P in the image field. The correction formulas are:     DX = (1-x-y + xy) DXa + (x-xy) DXb + (y-xy) DXc + (xy) DXD     DY = (1-x-y + xy) DYa + (x-xy) DYB + (y-xy) + DYC (xy) DYd where abc and d relate to the points ABC D. It has been shown that obtaining results DX and DY 12 bits could be achieved by calculating factors and correction DX DY eleven bits and calculating interpolation coefficients such as (1 - x - y + xy) 8-bit only. By simplifying the explanation we can say that getting fail tenth of a millimeter (12 bits for a 400 mm diameter image field) need to add at least two additional bits (must be as accurate as the quarter the finest gap to show). With the sign bit This leads to having to handle 15 bits.

However, in the invention usefully are measured only coordinate movements and not corrected coordinates. Considering that the maximum displacement would be met in the range of roughly 24 millimeters (which would be truly exceptional and would require actually a readjusting of the gamma camera) was able to realize that we could well gain a factor 400/24 ​​= 16, which is equivalent to four bits. It follows that it is no longer necessary to manipulate 15-4 = 11 bits (one sign bit, 8-bit dynamics and 2 fractional bits). And calculating interpolation coefficients on 8 bits only, because of a similar reasoning can be obtained results DX and DY 12 bits. Using to perform interpolation operation summing accumulators working on sixteen bits can easily implement the required accuracy. The latter explanation clearly shows that the accuracy is in fact obtained without increasing the computational power implementation with respect to what was described in the prior art.