DE102011115577A1 - Method for avoiding artifact in iterative reconstruction of object such as apple, involves determining start model of object function under back projection of detected gray values of projections for all prospects in object space - Google Patents

Abstract

The invention relates to a method for calculating a numerical model of the three-dimensional distribution of x-ray attenuation values for an X-rayed light from different directions, in particular a method for processing a plurality of electronically recorded X-ray images and for reconstructing the X-ray attenuation illuminated object, such as a living patient, by means of digital tomosynthesis.

Description

The invention relates to a method for calculating a numerical model of the three-dimensional distribution of x-ray attenuation values for an object illuminated with x-ray light from different directions. In particular, the invention relates to a method for processing a plurality of electronically recorded X-ray images and for reconstructing the transilluminated object, for example a living patient, by means of digital tomosynthesis.

wobei Ω⊂IR³ ein das durchleuchtete Objekt enthaltendes Raumgebiet (i. F. Objektraum) ist, die z-Achse entlang der Strahlrichtung verläuft und das eindimensionale Integral gerade die gesuchten – im Prinzip kontinuierlich verteilten – Abschwächungswerte f(z) des Objekts integriert. Jeder Pixel des Detektors erfasst die integrierte Abschwächung eines anderen Objektbereichs, d. h. einer anderen Strahllinie durch das Objekt, so dass die dreidimensionale Verteilung der Abschwächungswerte (i. F. bezeichnet als Objektfunktion f: Ω → IR) durch ihre Linienintegrale in den Grauwerten jedes Röntgenbildes zweidimensional codiert ist. Jedes Röntgenbild ist eine Projektion des Objektes. Stehen Projektionen für das gesamte Winkelintervall [0, π) in ausreichender Dichte zur Verfügung, so hat man die Radon-Transformierte der Objektfunktion aufgezeichnet. Aus ihr lässt sich die Objektfunktion selbst eindeutig rekonstruieren, was die Aufgabe der Computertomographie darstellt.Computed tomography (CT) is used in medical diagnostics, for example, when examining patient tissue (eg breast tissue) by means of X-ray light. The tissue is thereby irradiated along different directions, and the X-ray light is attenuated to different extents in the tissue types present along the respective beam direction (in particular by absorption and Compton scattering) and recorded in transmission on an electronic detector. The detector has a plurality of pixels arranged in a surface (also: point detectors) for spatially triggered measurement of the transmitted radiation intensity. The detector is hereinafter referred to as a flat detector and includes both flat-panel and curved pixel arrays. The radiation intensity I results from the irradiated intensity I ₀ according to

where Ω⊂IR ^{3 is} a spatial region containing the transilluminated object (i.f., object space), the z-axis is along the beam direction, and the one-dimensional integral just integrates the sought-after attenuation values f (z) of the object. Each pixel of the detector detects the integrated attenuation of another object area, ie another beam line through the object, so that the three-dimensional distribution of the attenuation values (i.f., referred to as object function f: Ω → IR) are two-dimensional in their gray line values of each x-ray image due to their line integrals is coded. Each x-ray image is a projection of the object. If projections are available for the entire angular interval [0, π) in sufficient density, then the radon transform of the object function has been recorded. From this, the object function itself can be clearly reconstructed, which is the task of computed tomography.

For various reasons (eg time savings, dose reduction, equipment restrictions) sometimes only a few projections of an object are available. For example, the object can only be transilluminated at very few angles or only at angles from a restricted angular range (in the first instance perspectives). Then the task of reconstructing the object function is an "ill-posed problem" that can not be solved unambiguously.

However, it can often determine approximate solutions that z. B. are sufficient for medical interpretations. Methods that calculate such approximate solutions from a limited selection of projections are summarized under the term "Digital Tomosynthesis (DT)".

• a. Ermitteln eines Startmodells der Objektfunktion unter Rückprojektion der detektierten Grauwerte der Projektionen für alle Perspektiven in den Objektraum,
• b. Numerisches Simulieren der Grauwerte wenigstens einer Projektion unter wenigstens einer Perspektive durch Vorwärtsprojektion des zuvor ermittelten Modells auf die wenigstens eine Detektorfläche,
• c. Vergleich der simulierten mit den detektierten Grauwerten und Ermitteln der Abweichungen,
• d. Aktualisieren des zuvor ermittelten Modells durch Addieren eines die Abweichungen korrigierenden Terms,
• e. Wiederholen der Schritte b bis d bis zur Konvergenz.

A subset of the known DT methods uses iterative algorithms. These iterative methods conceive the reconstruction problem as an optimization task and generate intermediate results for the object function sought, which are improved in each further iteration step until a convergence criterion is met. An iterative algorithm typically includes the steps:

• a. Determining a start model of the object function by back projection of the detected gray values of the projections for all perspectives into the object space,
• b. Numerically simulating the gray values of at least one projection under at least one perspective by forward projection of the previously determined model onto the at least one detector surface,
• c. Comparison of the simulated and the detected gray values and determination of the deviations,
• d. Updating the previously determined model by adding a deviation correcting term,
• e. Repeat steps b to d until convergence.

The gray levels of the projections are usually defined as a logarithm of the ratio of the radiated to the detected intensity, so that they are directly proportional to the integrals of equation (1). Their simulation from the model of the object function can therefore be done via a linear mapping.

For example, if i = 1, ..., N counts the pixels of the at least one area detector and k = 1, ..., P the available perspectives of the projections, then G _i ^{(k) should be} the measured gray value of the projection to the perspective k on the pixel i designate. All readings can be written to a column vector of length N + P, such as

The space region Ω is subdivided into voxels with indices j = 1, ..., M, so that the true object function can be understood as a step function on the voxels and represented by a column vector F of length M with the components f _j . Then, the forward projection mentioned in step b above is nothing more than a matrix multiplication of the vector F with a (N + P) × M system matrix R with constant coefficients resulting from the geometrical conditions of the transillumination of the object, inter alia from the arrangement of Radiation source and detector and the opening angle of the radiation result. This two-dimensional system matrix with N + P rows and M columns is known.

The task of digital tomosynthesis is to determine the unknown components of F , although the system of equations (3) is underdetermined and there are infinitely many solutions. Each of these solutions is a possible model of the object.

und können ebenfalls in das Objektvolumen zurückprojiziert werden, um ein Update des Modells etwa durch Subtraktion der zurückprojizierten Abweichung zu bestimmen. Der n-te Iterationsschritt lautet dann F ⁽ⁿ⁾ = F ^(n-1) – λR ^T δ ^(n-1), (6) wobei δ ^(n-1) wie in Gleichung (5) stets mit F ^(n-1) verknüpft ist. Der reellwertige Parameter λ wird als Relaxationsparameter bezeichnet und vom Anwender des Algorithmus aus dem Intervall 0 ≤ λ ≤ 1 (oder 2 in manchen Publikationen) gewählt. Der die Abweichung korrigierende Term aus Schritt d. oben wird üblich mit dem Relaxationsparameter multipliziert, wie in Gleichung (6) angegeben.In order to iteratively approximate a solution, the measured values G are first projected back into the volume Ω (more precisely, to the voxels located there) by multiplication with the transposed system matrix. The result is a starting model of the object function, such as F ⁽⁰⁾ = R ^T G , (4) However, when inserted into equation (3) - ie in the new forward projection - the vector G only approximately reproduced. The deviations are given by

and may also be back-projected into the object volume to determine an update of the model by, for example, subtracting the back-projected deviation. The nth iteration step is then F ⁽ⁿ⁾ = F ^(n-1) -λ R ^T δ ^(n-1) , (6) where δ ^{(n-1) is} always associated with F ^(n-1) as in equation (5). The real-valued parameter λ is called a relaxation parameter and selected by the user of the algorithm from the interval 0 ≦ λ ≦ 1 (or 2 in some publications). The deviation correcting term from step d. above is commonly multiplied by the relaxation parameter as given in equation (6).

As is known, the naive implementation of the iteration method outlined so far is not recommended because the required computing time and the memory requirement are considerable. The very large system matrices are only very weakly occupied, so that there is much room for numerical optimizations. In particular, it is usually considered convenient to use the different perspectives of the projections separately during the iteration, usually in succession.

To further accelerate and / or improve the iterative reconstruction, the prior art knows, for example, the Fourier transform of the measured values and system matrices and the subsequent iteration in the Fourier space (FIG. EP 1 793 347 A2 ), the introduction of voxel-dependent weight functions adapted in each iteration step ( US 2008/205729 A1 ), filtering and / or smoothing of measured values and / or highlighting similar ranges of measurements in the individual projections when creating the start model of the iteration ( US 6,707,878 B2 and US 6,973,157 B2 ) or the selective iteration of only one "region of interest ", which is to be determined by the user after a first, quickly calculated representation of the start model for the more intensive - intensive evaluation ( US 2008/187090 A1 ).

In fact, the expressiveness of the starting model is generally high, so that one can sometimes be satisfied with its computation from the projection data - as it were with the 0th iteration - as a result of the reconstruction. Among the optimization methods that already start with the reconstruction of the start model (and which can also be repeated if necessary in subsequent iteration steps), two are particularly noteworthy.

The US 2003/072478 A1 discloses that it may be expedient, after the backprojection of all measured data (vector G ) on the voxels of the object space Ω not immediately calculate the starting model according to equation (4), but rather cancel the matrix multiplication in (4) before the summation on the perspectives and instead submit the contributions of all perspectives to each voxel as an ensemble to a nonlinear operator. For example, this may be the minimum operator, ie each voxel is assigned the smallest contribution of all perspectives that can be assigned to it during the backprojection.

The simple idea behind it is that even if back projections of almost all perspectives assign a high absorption value to one voxel, but one does not, the latter is probably correct, and the absorbing structure of the object that causes the measurements will be more likely to be in another, probably in the To be found nearby voxels. The presence of an absorbent structure is indicated by the "unanimity of the vote", ie by a high contribution of all perspectives to a voxel.

The US 2003/072478 A1 but also provides other possibilities for the non-linear operator, such as maximum, median or filtering for a given threshold. The totality of the values assigned to the voxels by the non-linear operator is used either as a result of the reconstruction or possibly as the starting model of an iteration.

The US 2009/060310 A1 similarly, it argues that for each voxel of the model, it requires the consistency of the backprojection of the different perspectives ("consistency requirement of the backprojection process", para. 0006) and causes it to interfere with the data. The intervention takes place after the back-projection of the measured gray values into the object space, by first adding to each voxel j (j = 1, ..., M) the contributions of the perspectives B ^(j) _k (k = 1, ..., P). after which a sorting and evaluation of the ensemble B ^{(j) is performed} for each fixed j and then such ensemble values B ^(j) _k satisfying the given criteria of significance (A1) are identified and replaced by other values , The significance conditions are explained in such a way that they determine and limit the deviation of individual ensemble values from the majority of the (possibly remaining) ensemble B ^(j) , for example the distance to the mean value or the position relative to the distribution width. In short, the method aims to detect outsiders in the ensembles B ^(j) and to exclude them from the further calculation for the reconstruction of the object model.

The invention now has the task of proposing an improvement of the iterative methods for digital tomosynthesis, in particular for better artifact avoidance.

• a. Ermitteln eines Startmodells der Objektfunktion unter Rückprojektion der detektierten Grauwerte der Projektionen für alle Perspektiven in den Objektraum,
• b. Numerisches Simulieren der Grauwerte wenigstens einer Projektion unter wenigstens einer Perspektive durch Vorwärtsprojektion des zuvor ermittelten Modells auf die wenigstens eine Detektorfläche,
• c. Vergleich der simulierten mit den detektierten Grauwerten und Ermitteln der Abweichungen,
• d. Aktualisieren des zuvor ermittelten Modells durch Addieren eines die Abweichungen korrigierenden Terms multipliziert mit einem Relaxationsparameter,
• e. Wiederholen der Schritte b bis d bis zur Konvergenz,

dadurch gekennzeichnet, dass nach der Rückprojektion der detektierten Grauwerte der Projektionen für alle Perspektiven in den Objektraum für jeden Voxel wenigstens eine die Unähnlichkeit des Ensembles der in den Voxel für alle Perspektiven zurückprojizierten Grauwerte indizierende Maßzahl berechnet wird und weiterhin jeder Maßzahl durch eine stetige, streng monoton fallende Funktion ein Wert zugeordnet wird, mit dem der Relaxationsparameter bei der Durchführung der Aktualisierung in Schritt d. wenigstens für jeden der Voxel separat initialisiert wird. Die Unteransprüche geben vorteilhafte Ausgestaltungen an.The object is achieved by a method for artifact avoidance in the iterative reconstruction of an object illuminated by X-ray object in an object space, wherein an object function is determined on a the object space dividing set of voxels from at least one flat X-ray detector under different perspectives measured gray value images of the projections of the object , with the steps:

• a. Determining a start model of the object function by back projection of the detected gray values of the projections for all perspectives into the object space,
• b. Numerically simulating the gray values of at least one projection under at least one perspective by forward projection of the previously determined model onto the at least one detector surface,
• c. Comparison of the simulated and the detected gray values and determination of the deviations,
• d. Updating the previously determined model by adding a deviation correcting term multiplied by a relaxation parameter,
• e. Repeating steps b through d until convergence,

characterized in that after the back projection of the detected gray values of the projections for all perspectives into the object space for each voxel at least one of the dissimilarity of the ensemble of the in the Voxel is calculated for all perspectives backprojected gray scale indicative index and further each measure is assigned by a continuous, strictly monotonically decreasing function a value with which the relaxation parameter when performing the update in step d. at least for each of the voxels is initialized separately. The dependent claims indicate advantageous embodiments.

The steps and calculations of the above method and of the method defined in the claims can be carried out-advantageously automatically or semi-automatically-preferably completely or at least partially with the aid of a suitably adapted or programmed data processing device.

The effect of the invention is illustrated below with reference to example images. Show

1a) A photograph of an apple prepared with metal needles, wherein the needles are arranged substantially in a plane which cuts the apple in the middle; b) sectional image (in needle plane) of the result of a DT reconstruction from X-ray images of the apple according to the prior art;

2 Sectional images of the apple analogous to 1b) reconstructed according to the method of the invention;

3 DT reconstruction of a human hand (selected section plane) according to the prior art;

4 DT reconstructions of the hand 3 from the same X-ray images according to the method of the invention.

The invention is based on US 2009/060310 A1 as the closest prior art and makes particular the idea that the contributions of the P perspectives in the form of ensembles B ⁽ in the voxelwise reconstruction of the searched object model after the backprojection - at least for the start model, but preferably also in later iteration steps) ^{j) should be} noted for each voxel j and examined for the similarity of their components.

In contrast to US 2009/060310 A1 should not be intervened in the data.

Rather, it is provided according to the invention to modify the iteration process by setting a separate relaxation parameter λ = λ _j for at least one voxel. The idea is that voxels with similar values can be iterated relatively robustly in ensembles B ^(j) - that is, with high agreement of all perspectives - whereas voxels with dissimilar values in ensembles B ^(j) tend to iterate with the formation of artifacts must be expected in the model.

It should be emphasized that the relaxation parameter has its origin in numerical mathematics (here: Kaczmarz method for the approximate solution of systems of equations), but not in tomography. Although one knows the possibility for the possible adjustment in the course of an iteration, but there are hardly any rules for its initialization. However, it is known that the initialization can definitely influence the iteration result.

The inventors have recognized that for use in iterative DT, the initial tomographic conditions should already be coded in the relaxation parameter in order to bring about a reduction of artifacts during the reconstruction. For this purpose, the relaxation parameter is to be expanded at least to all voxels of the object model, i. H. as a vector or array with voxel index j, and also to initialize at least for each voxel individually.

The iteration of the object model then takes place according to the invention in analogy to equation (6) on the individual voxels in component notation as follows:

If one writes the deviations analogously to equation (5), the result is:

The ensemble values B ^(j) _k are the detected gray levels projected back into the object space, and the ensemble values C ^(j) _k are also backprojected forward projections of the intermediate values of the iterated object model for voxel j and perspective k as defined by the parenthesized expressions in Equation (). 8th). No assumption is made about the concrete shape of the system matrix R or R ^T. The prior art knows several variants of how the projections can be defined. As long as forward and backward projection are writable as linear mappings, equation (8) represents the most general formulation of the iteration rule.

The invention proposes to make the parameters λ _j at least dependent on the B ^(j) _k . As a general rule, each λ _j for fixed j should be chosen so that the smaller the dissimilarity of B ^(j) _k for the P perspectives, the smaller it is. Preferably, a measure of the dissimilarity of B ^(j) _{k is} introduced, assuming on each voxel a value between zero and one, where zero indicates the equality of B ^(j) _k for all k = 1, ..., P , It is possible that the value one will never be accepted (see below).

The dissimilarity of the ensemble values (abbreviated as D in the following) is mapped according to the invention via a continuous, strictly monotonically decreasing function ψ, and finally one identifies λ _j = ψ (D ^(j) ), (9) where D ^{( j)} denotes a measure of the dissimilarity of the ensemble B ^(j) . Preferably ψ is a function with the properties ψ: [0, 1] → [0, 1] with ψ (0) = 1 and ψ (1) = 0.

In principle, there are as many possibilities as possible to define the dissimilarity. Common to all expedient definitions, however, is that for given j they comprise the determination of a value Ξ _j representing the ensemble B ^(j) and the determination of a distance measure of all B ^(j) _k to Ξ _j . Preferably, Ξ _{j is identified} with one of the ensemble values B ^(j) _k , such as by selecting the smallest or largest ensemble value or the median of the ensemble

Alternatively, Ξ _{j can} also be determined as the arithmetic mean of B ^(j) _k .

According to the invention, the calculation of the distance measure comprises the formation of a norm of the differences of B ^(j) _k to Ξ _j , the Euclidean (L2) standard being used by way of example below and being noted by simple absolute value strokes.

und erhält positive Maßzahlen aus dem Intervall [0, 1], wenn min(B^(j) _k) ≤ Ξ_j ≤ max(B^(j) _k) gilt. Man kann dann bevorzugt den arithmetischen Mittelwert bilden zur Bestimmung von

der folglich auch aus dem Intervall [0, 1] stammt, wobei der Wert Eins offenbar nie angenommen werden kann.For example, you calculate

and obtains positive measures from the interval [0, 1] when min (B ^(j) _k ) ≤ Ξ _j ≤ max (B ^(j) _k ). One can then preferably form the arithmetic mean for the determination of

which therefore also comes from the interval [0, 1], whereby the value one can obviously never be assumed.

The calculation of the voxel-dependent relaxation parameters according to equation (9) uses the results of equation (12). The strictly monotonically decreasing function ψ is in the simplest case as ψ (x) = 1 - x (13) established. Sometimes, however, it is advantageous if ψ falls faster than linear, so that, for example, too ψ (x) = 1 - x / 1 + x (14) comes into consideration. It has even proved to be particularly advantageous to introduce a throttle factor α, which enforces an even greater decrease of ψ for α> 1, and to make this throttling factor inter alia also dependent on the distribution of d ^(j) _k from equation (12). ψ (x) = (1 - x / 1 + αx) ^α (15)

A steeper curve of ψ leads to smaller relaxation parameters with greater dissimilarity, ie. H. to a more "careful" iteration, especially on those voxels that can give rise to artifact formation. On the one hand, the choice of α takes into account the degree of dissimilarity on how many voxels is to be determined and, on the other hand, the mean dynamic range present in the measured data must be taken into account.

The average contrast can be determined, for example, by the fact that the totality of all measured values B ^(j) _k projected back into the object space Ω is first ordered by size for each voxel j, and hereafter an average value <B _max > for the largest and <B _min > for the smallest values is determined as the arithmetic mean over the voxels.

The difference between these values can be considered as medium contrast. If the set of measured values is high in contrast, ie the calculated average contrast is large, then voxels with high dissimilarity will as a rule also have a high degree of dissimilarity D ^(j) , and the steepness of the function ψ, and consequently the throttling factor α, should be reduced.

und ein erfindungsgemäß zweckmäßiger Drosselfaktor ergibt sich gemäß

zur Verwendung in Gleichung (15). Mit β ist ein Proportionalitätsfaktor bezeichnet, den der Nutzer nach Bedarf wählen kann. Die übliche Wahl ist β = 1.In contrast, the throttle factor is to be increased if the data on many voxels show relatively large dissimilarity measures. It has proved to be advantageous to determine the centroid of mass of the dissimilarity measures according to equation (11) and to use it to establish α. For this purpose, preferably all are sorted by size and assigned according to their size discrete intervals and each counter a counter, which is increased by one, as soon as an assignment in this interval takes place. In other words, a histogram of d ^(j) _{k is} created, which assigns the frequencies (here: interval centers d _i , i = 1,..., V with V as the number of intervals) to frequencies n _{i of} the occurrence. The numbers of the measure allocations are noted for all intervals; a size distribution of the measures is determined. The focus of the distribution is then

and an appropriate throttle factor according to the invention results according to

for use in equation (15). Β denotes a proportionality factor that the user can choose as needed. The usual choice is β = 1.

The method described so far produces a voxel-dependent initialization of the relaxation parameters λ _j for the iterative reconstruction of the object function in digital tomosynthesis. The method is parameter-free and is based exclusively on the measured values which are examined immediately after the first backprojection into the object space. In particular, the dissimilarity of backprojected measured values from different perspectives is quantified for the individual voxels of the object model and thus forms the basis for the determination of the λ _j .

It can be assumed that other iterative DT optimization methods that do not address the relaxation parameter can be easily combined with the method described herein. Negative effects of such combinations have not been recognized.

The method described here can be used repeatedly even during the iteration. In particular, the determination of dissimilarity may also be made for ensembles of backprojected forward projections of the intermediate values of the iterated object model, C ^(j) _k, in equation (8). The relaxation parameters can then be calculated exactly as described and used in the following iteration step.

Furthermore, it may be advantageous to set a separate set of relaxation parameters for each individual perspective. This is an option if the iteration is designed in such a way that the object model should be updated in each iteration step with only one perspective (ie one of the P x-ray images). Suitable relaxation parameters are obtained by λ _k ^(j) = ψ (d _k ^(j) ), (18) which takes the place of equation (9). The statements of equations (11) and (13) to (17) remain unaffected.

Finally, the improved artifact avoidance in DT reconstruction according to the method of the invention will be demonstrated by way of examples.

In the first example, an apple is spiked with metal needles by hand so that the needles penetrate approximately centrally into the apple and all lie approximately in one plane (needle level), as in 1a) shown. X-ray images of the thus prepared apple are recorded from different perspectives, with the X-ray source being moved from left to right in all the images shown.

From the gray values of the X-ray images, an object model of the apple is calculated using an iterative standard method (here: Simultaneous Algebraic Reconstruction Technique, SART). The starting model is initialized to zero, there is an update in each iteration step with the data of exactly one perspective in random order, and each perspective is used only once.

The voxelization of the object space is such that 60 slices, each 1 millimeter thick, are modeled parallel to the (anticipated) needle plane. In 1b) is the result of the prior art reconstruction in level 29. Some needles - z. B. highlighted by ellipses in the picture - are exactly in the needle level, while others are clearly assigned to another level and some apparently also have an angle with respect to the needle plane. The image shows typical DT artefacts that are indicated by arrows: Shadow artifacts in the lower needle environment, ghost copies of the upper needles, needle-shaped shadows on the side needles, and finally quite pronounced contour artifacts around the entire apple.

The reconstruction of the same X-ray images with SART and the method according to the invention leads to the images in 2 , Also shown here are the calculated planes 29, and the markings (arrows, ellipses) are for easier comparison 1b) Have been transferred. For the calculation of 2 Relaxation parameters are initialized for each voxel and for each perspective according to equations (11) and (18), where 2a) the linear function according to equation (13) and in 2 B) the convex function of equation (15) with α = 2 is used.

It is by comparing the 2 With 1b) readily apparent that the DT artifacts are significantly reduced and z. T. are virtually eliminated. The use of the convex function of Eq. (15) seems to give somewhat better results if metallic objects are the cause of the artifacts.

The second example is the reconstruction of a human hand, for which an object model of 44 levels according to the state of the art (also SART) is first created. A plane of this model, in which the proximal phalanges of the middle three fingers, the metacarpal bone of the little finger and the carpal bones, are focused (marked by ellipses), is in 3 to see. Again, arrows point to areas with pronounced artifacts.

4 Finally, the results of the reconstruction according to the invention show markedly reduced artifacts. Again, the relaxation parameters are initialized for each voxel and each perspective, where for 4a) the function according to equation (13) and for 4b) which is used according to equation (14). Both provide very similar results here.

The above steps and calculations can be carried out advantageously - preferably automatically or semi-automatically - completely or at least partially with the aid of a suitably adapted or programmed data processing device.

QUOTES INCLUDE IN THE DESCRIPTION

This list of the documents listed by the applicant has been generated automatically and is included solely for the better information of the reader. The list is not part of the German patent or utility model application. The DPMA assumes no liability for any errors or omissions.

Cited patent literature

• EP 1793347 A2 [0012]
• US 2008/205729 A1 [0012]
• US 6,707,878 B2 [0012]
• US 6973157 B2 [0012]
• US 2008/187090 A1 [0012]
• US 2003/072478 A1 [0014, 0016]
• US 2009/060310 A1 [0017, 0026, 0027]

Claims (10)

A method of artifact avoidance in the iterative reconstruction of an X-ray illuminated object in an object space, wherein an object function is determined on an object space dividing set of voxels from gray level images of the object's projections measured with at least one X-ray area detector under different perspectives, comprising the steps of: a , Determining a start model of the object function by backprojecting the detected gray values of the projections for all perspectives into the object space, b. Numerically simulating the gray values of at least one projection under at least one perspective by forward projection of the previously determined model onto the at least one detector surface, c. Comparison of the simulated with the detected gray values and determination of the deviations, d. Updating the previously determined model by adding a deviation correcting term multiplied by a relaxation parameter, e. Repeating steps b to d until convergence, characterized in that after the backprojection of the detected gray values of the projections for all perspectives into the object space for each voxel, at least one measure indicating the dissimilarity of the ensemble of the gray values projected back into the voxel for all perspectives is determined and further assigning to each measure a value by a continuous, strictly monotonically decreasing function, with which the relaxation parameter in performing the update in step d. at least for each of the voxels is initialized separately. Method according to claim 1, characterized in that the at least one measure calculated for each voxel is determined from the interval [0, 1] and the continuous, strictly monotonically decreasing function assigns each measure from the interval [0, 1] to the interval [0 , 1] depicts. A method according to claim 2, characterized in that the continuous, strictly monotonically decreasing function by ψ (x) = 1 - x is given. A method according to claim 2, characterized in that the continuous, strictly monotonically decreasing function by ψ (x) = (1-x / 1 + αx) ^α given with a throttle factor α> 1. Method according to one of the preceding claims, characterized in that the measure for the dissimilarity of the ensemble of the gray values projected back into a voxel is calculated as a norm of the differences of the ensemble values representing the ensemble value, divided by the largest ensemble value and averaged over the ensemble. Method according to one of claims 1 to 4, characterized in that the steps b. to d. in a predetermined iteration step with a predetermined projection to a predetermined perspective, and that for each voxel and for each perspective a difference in the ensemble of the gray values indicative of the voxel for all perspectives is calculated, and furthermore each measure is represented by a continuous, strict monotonically decreasing function is assigned a value with which the relaxation parameter in the implementation of the update in step d. is initialized separately in the predetermined iteration step for each of the voxels and for the predetermined perspective. Method according to one of claims 1 to 4 or 6, characterized in that the for a predetermined perspective, the dissimilarity of the ensemble of gray scale indicative of the voxel backprojected for all perspectives as the divided by the largest ensemble value norm of the difference of the predetermined perspective assigning Ensemble value is calculated to a value representing the ensemble. Method according to one of claims 5 or 7, characterized in that the value representing the ensemble is selected from one of the following options: the maximum, the minimum, the median or the arithmetic mean of the ensemble values. A method according to claim 7, characterized in that the measures indicating the dissimilarity of the ensemble are assigned to predetermined intervals for all perspectives and for all voxels according to their size, each measure being assigned exactly once to one of the intervals and the numbers of measure assignments for all intervals is determined as the size distribution of the measures. A method according to claims 4 and 9, characterized in that the throttle factor is proportional to the ratio of the centroid of the size distribution of the measure indicating the dissimilarity of the ensemble and the mean contrast of the gray scale images.

Images