FI116324B - Method and arrangement for producing a tomography image - Google Patents

Description

1U524

The invention relates to the production of tomographic images of irradiated objects in imaging systems in which the transmitted radiation of the object is detected by a digital detector having a specific spatial resolution. In particular, the invention relates to utilizing statistical inversion in image production.

In general, tomography imaging is the imaging of substantially internal structures of a solid object by observing and recording differences in the behavior of certain transmitted waves as they strike those structures. One typical application is medical x-ray tomography, where the subject is a live organism or part thereof, the radiation waves used are X-ray paths ranging from a few keV to tens or even about 100 keV and the purpose of the imaging process is to make diagnostic observations are not visible from above. Other tomography applications include various industrial processes where it is useful to obtain information about the intrinsic properties of a given raw material or product. For example, log tome-20 graphics are designed to examine logs prior to sawing so that each log can be sawn into planks in an optimal manner.

. ; Figure 1 illustrates the basic principle of two-dimensional low-angle computed tomography.

'; The radiation source, shown in the figure at two exemplary locations 102 • 25 and 103, irradiates the object 101 from a limited number of different directions. A spatially sensitive linear detector shown at two positions 104 and 105, respectively; . 'Forms a damping profile for each exposure. Figure 1 assumes that the quoted areas within the object 101 attenuate radiation to a greater degree than the target effective mass. The damping profiles and corresponding irradiation angles are fed to a computer ··· 30 machine 106, which calculates, by a mathematical method, how the damping elements are to be located inside the object 101 to produce the damping profiles. The result of the calculation that substantially corresponds to the two-dimensional intersection of • object 101; 1 damping coefficient map, shown on screen 107.

: ': 35 A number of different mathematical reconstruction methods have been proposed to convert the obtained measurement results into an object image. Known and successfully used reconstruction algorithms include filtered back projection (FBP), least squares method, maximum likelihood estimation, algebraic reconstruction (ART), simultaneous iterative ART (SIRT), and multiplicative ART (MART). All of these work reasonably well if there is a large number of views from different directions and the views cover a wide area around the subject. However, multiple exposures involve exposing the subject to high doses of ionizing electromagnetic radiation, which, in the case of medical x-ray tomography, is contrary to the general goal of keeping the dose of radiation received at a low level. In many cases, it is physically impossible to freely choose the position of the subject and the imaging arrangement relative to one another so as to achieve a large number of different angles. The construction of large rotary arrangements 10 or the use of multiple parallel measuring sensors often makes the arrangement more complex and prone to failure. It also increases the cost of manufacturing.

A promising advantage over the conventional methods described above has been obtained by utilizing 15 statistical inversions. As a mathematical method, it is not new, but has long been considered too demanding for computational resources to be applicable in practical applications. In particular, the present invention relates to a so-called Bayesian inversion based on the following formula: 20 p (4,) - w (1) where D {x \ m) is a probability density function corresponding to image data, Dpr (x) is so; called a priori, that is, a mathematical model of certain properties of an object, D (m \ x) is a * ;, / absorption pattern that associates a vector of the object a with the projection result vector ra, ··· '25 and D (m) is a normalization constant. The significant features of Bayesian inversion are: first, that image representation a is not an unknown deterministic

* <I

constant, as in normal maximum likelihood estimation, but a random variable, and second, that a priori Dpr (x) is independent of the measurement result m.

· **, 30 The constant D (m) does not depend on the image, which means that it does not have to be taken into account when searching for the most suitable image data estimate. The latter task can be written: 'in the form:' j '. a = arg max [z] (ra | x) D (a)] (2) ·; *: 35 116324 3 which means that the result given as the final image is the estimate x that maximizes the right expression of equation (2). The purpose of the image reconstruction process is thus to find the right equilibrium which, on the one hand, provides a good correspondence to the measured data p (m \ x) and on the other hand does not overly violate the a priori assumptions contained in a priori Dpr (x). The solution of the Bayesian inversion problem involves the creation of random values from the so-called posterior distribution.

Priori is a set of rules about what can and cannot appear in an image. If these rules are too strict, it is possible that essential information will be lost. 10 Ideally, if an absolutely true image of the subject could be obtained, it would not be "fined" a priori at maximizing the expression in equation (2). Classical examples of priests can be found in the so-called. from the Gaussian or normal priority class. However, Gaussian priorities are not the best choice for medical X-ray imaging or certain other important applications of computed tomography, because they assume that the absorption characteristics of the subject change steadily from one part of the subject to another. For example, in parts of the human body, sharp transitions occur between bone and soft tissue, whereas, for example, within a constant density bone, the absorption properties remain substantially unchanged. Other major applications, such as geological sample analysis, also exhibit similar structural features. The use of a Gaussian prior would lead to degradation of the image quality by blurring the edges between the two absorbent regions of the subject and,. ·. creating variation in areas where it should not occur. Adjusting the number of views li- 2 usually reduces the adverse effect of the Gaussian priorities, but

I t I

its use is not always recommended or even possible for other reasons.

': .2 25

Total Variation Priorities are an example of priorities that are much better suited for medical imaging and other applications where there is a steep displacement of the ... ... · sea buckthorn between different absorption areas. The total variation priority is a function whose exponent contains pixel difference in pixels of adjacent pixels; i · 30 yen. The total variation priorities tend to steer the computation toward solutions of particular regions such that the absorption properties of each region are substantially constant, while at the boundary of regions there may be a distinct discontinuity. The disadvantage of total variation priorities is that, mathematically, they ''; · 'Are more complex to handle than, for example, Gaussian priorities, which slows down the image reconstruction process.

i * * »{» t

There are other problems with known methods of applying Bayesian inversion to tomographic reconstruction. The basic principles of all digital image processing include the use of a grid, whereby in a two-dimensional case the image is made up of a plurality of pixels. A three-dimensional object is modeled on a three-dimensional lattice consisting of volumetric pixels or voxels. Each pixel or voxel is associated with a value, for example, in X-ray tomography, the value of the absorption coefficient of 5, which at the beginning of the calculation is an unknown variable. Together, unknown variables form unknown vectors and matrices. The recon-structuring process thus involves a large number of matrix computations in which the required number of operations is a square or a cube of the number of pixels. Thus, it is difficult to obtain meaningful results within a reasonable time and with reasonable computation power.

It is an object of the present invention to provide a method and arrangement for applying statistical inversion to tomographic reconstruction, significantly reducing the problems associated with the prior art. In particular, it is an object of the invention to reduce the number of CPU operations performed in image reconstruction and thus enable the use of statistical inversion computed tomography with computing resources available in any treatment or industrial unit. It is a further object of the invention to maintain the basic versatility of statistical inversion so that the method and arrangement 20 of the invention can be used in a variety of tomography applications.

. The objects of the invention are achieved by the use of a priori consisting of parts which, in themselves, are easy to handle in computation. Other objects of the invention are achieved

>! I

'! utilizing the so-called. the rarity of the theory matrix to avoid unnecessary calculations.

s *:. * The process of the invention is characterized by what is presented in the independent. 'Method requirement in the character part.

The invention also relates to an arrangement characterized by what is stated in the characterizing part of the independent claim directed to the arrangement.

The invention further relates to a computer program product. The computer program product '· ·' according to the invention is characterized by what is represented by a computer program product; ': In the characterizing part of 35 independent claims relating to a product.

The first aspect of the invention is based on an insight into the properties of priori and their effect on computing resource requirements. One feature known per se, 5 1 16324, is that the Gaussian priorities also make the posterior distribution a Gaussian or normal distribution. This has the advantage that generating random values to be used in the image reconstruction process is relatively easy because good, fast and robust general-purpose algorithms are readily available to generate random numbers of some kind of normal distribution. Instead, the use of total variation priorities implies the need to perform a large number of numerical integrations during the image reconstruction process to form the poster theorems from which random values should be constructed. Performing integrations requires a very large number of processor operations, which in turn requires either supercomputer capacity or a very long time or both. The computational task is by no means facilitated by the fact that the coordinate region in which the integration is to be performed is necessarily very large, since it is necessary to have certainty that the probability density outside that region is essentially zero.

The first aspect of the invention involves treating the posterior distribution as fractions which themselves follow a distribution that makes it easy to construct random variables. As a prime example, the posterior distribution may be treated as a piecewise Gaussian distribution, but the invention does not exclude the use of other distributions that satisfy the same condition, namely the relatively easy generation of random numbers. When the posterior distribution is made up of such parts, random values are generated by drawing the relevant part. When a part is selected, the properties of that distribution are used to construct the value of the random variable.

► i · ·· *: Another aspect of the invention relates to the so-called. processing of the theoretical matrix. The theory matrix emerges in the exponential posterior distribution associated with the actual measurement, not a priori. Traditionally, prior art algorithms do not use the theory matrix itself, but a certain other matrix obtained by multiplying the theory matrix by its transpositions. This results in a situation where each iteration of the image reconstruction algorithm requires matrix times vector multiplication, which may contain 30 tens of millions of nonzero elements. Not only is the number of processor operations needed to perform multiplication directly, · but huge is also the * P memory space required for all the elements in those matrices and vectors.

i i t

According to another aspect of the invention, the theory matrix is not multiplied by its transpositions,. ·:, 35 but is used as is. This allows us to take advantage of the rarity of the theory matrix, which is a direct consequence of the so-called basic assumption used in the tomography system. teräväkeilamallista. In addition, another aspect of the invention is based on the fact that when a vector of pixel or voxel values is multiplied by a theoretical matrix, it is not necessary to perform a multiple matrix times vector multiplication. Suffice it to update the matrix times vector input for that coordinate. This avoids a significant amount of virtually unnecessary calculations.

The novel features considered to be characteristic of the invention are set forth in detail in the appended claims. However, the invention itself, its structure, its working principle, and its additional objects and advantages will be described below with reference to some embodiments and with reference to the accompanying drawings.

Figure 1 illustrates the basic principle of computed tomography, Figure 2 illustrates an example of a piecewise Gaussian conditional distribution, Figure 3 illustrates the relationships between illumination geometry and the theoretical matrix, Figure 4 illustrates the high level method steps of the invention, Figure 5 further illustrates the method 7, a schematic view of the device according to the invention and FIG. 8 a computer program product according to the invention.

Exemplary embodiments of the invention disclosed in this patent application are not to be construed as limiting the applicability of the appended claims. The term "comprising" is used in this patent application as an open-ended limitation. which does not exclude features not mentioned here. The features disclosed in the dependent claims may be freely combined with one another; unless otherwise stated.

25

As a background to the present invention, some known features of statistical inversion are briefly set forth below.

Suppose there is a so-called random variable with respect to x. posterioryakauma f (x). In order to find the most acceptable value of x, the Markov chain is formed by starting with »·» *, ··. of x, and use the general rule that each x; + 1 x, • 'is formed by taking it from the probability distribution d (t + iK) = p (t-T + i) · (3)

>> I I I I I

35 116324 7

Assuming now that / M =?) <** »(4) 5 then the ergodic assumption implies that tv (= 1 for any measurable function g and the number of elements of chain (jt;) N. 10 In other words, chain (* < ) can be used as a simulated sample of the posterior distribution fix) for example to calculate any of the higher moments of xn mean E fix), posterior variance Et (x1) - Et (xf or xn).

Two commonly used methods for forming the Markov chain, known as MCMCs (Markov Chain - Monte Carlo), are the so-called Metropolis-Hastings or Gibbs sampling. The MH method generates the entire new sample directly, while the Gibbs sample generates one new coordinate at a time. The MH method is relatively fast with all sorts of distributions, but it does not work with a large vector jc dimension. Gibbs sampling instead of 20 operates on arbitrarily large dimensions, but it has the disadvantage that it f t. T is fast enough in practical applications with only a few distributions. For additional constraints, such as the positivity requirement, Gibbs sampling is easy, • act; ·; like the MH method.

»* • ..: 25 Let's take a closer look at Gibbs sampling using the still-assumed posterior distribution j fix) but this time using the / V-dimensional vector variable random variable x, ,, _: whose elements are (jr !, x2, jc3, ... xfi ). This time, Markov's chain formation for converting xn to another vector y is denoted as:, ··, 30 1) Construct yj from the distribution /(χ^χ.,ΐΦΐ):: · 'D (yx) ~ f {yx, x2, x3, ..., ½) (6) •> <(7) then form y2 such that 35 116324 8 3) and continue until (8) 5 This procedure defines a random transformation x -> y. The transition density function can be denoted p (x, y) = D (y | jc), (9) 10 and it is easy to see that the criterion of equation (4) is satisfied. Thus, Gibbs sampling essentially generates random numbers from a one-dimensional distribution. With certain distributions, such as even, normal, exponential and some other distributions, this is easily accomplished. In the general case, however, the independently known fact that the density of the function F1 (r) is / if r is a uniformly distributed random number and F is a distribution function fin must be used. When performing a Gibbs sampling with a general posterior distribution /, it is necessary to compute F, that is, integrate the function f over such a large area that the probability density outside of it is certainly substantially zero. This kind of numerical integration can be extremely time consuming with the best computers and algorithms if the problem is large.

20

Let us now look at the "easy" distributions, especially the normal distributions. . Initially, it is assumed that the posterior distribution / corresponds to a given normal distribution with; '[, mean xq and variance Q 1:

* I I

25 f (x) ~ N (x0, e- '). (10) * *: 'Denote • · Y = {x2, ..., xJ (11) 30. Q \ y (12) * »: '35 Then the conditional density of x ^ is given by:

! I

(13) 116324 9

If the posterior density / is obtained by inverting the linear direct theory ιη = Αχ + ε, (14) 5 where m represents the observation results, A is a theoretical matrix, λ: is a vector of unknown random variables and εοη is an error term, then density Dix) is £) (. *) = exp - - (m- Αχ) τΣ ~] (ηι- Αχ) 1 w (15) = exp - (xTQx-2xTv + C)

V 2 J

j j Tl 10 where Q = A Σ A, v = A Σ 'm, Σ is the error covariance matrix (which is diagonal, assuming that the measurements are independent in this sense) and C is not dependent on x. Thus, the centers of one-dimensional conditional distributions can be calculated without calculating the center of the entire posterior density. Calculating the latter would be time consuming if the number of unknowns is large.

15

Finally, some general forms of known priorities are presented. The Gaussian (atypical) prior is of the form f (χ-χ) 2); Dpr (x) ~ exp - £ '2J, (16) /' _ ((ij ^ N · "ij j: 20 '·': where V is a set of indices (i, j) such that x (and Xj are adjacent points and '.' djj is a regularization error, which here corresponds to the standard deviation and controls the relationship between measurement and re-; \: gularization. Providing indexes for a regularization error is only possible if its variation between lattice points is allowed; In equation (16), the regularization error is treated as a constant parameter. j · The fact that the priori is "atypical" means that the integral of the priori itself has a peripheral value. '' ·.

; ' (aij J

: 30 'The total variation priority can have a positive real number 5 as an exponential exponent of the difference, ie | *, · -x: y |' in summable terms, and not just | jc ,. - x ^.

116324 ίο

A first aspect of the invention relates to the way in which the use of a total variation type priority can be combined with a very fast and efficient generation of random variable values. Assuming a total variance priority and a particular fixed vector * o = (* io. * 20. ···). it can be shown that the conditional distribution 1-½.-½ ....) has 5 forms / \ f if 2 f * Y) | λ:, - λ:; · 0Π F20 '-½)> · · ·) ~ exP - TX \ Q \\ ~ ^ -X \ Vi - ^ ~ (18) \ V / = 2 JJ (1 J ^ N- a \ j Here, Qu and Q12 are the actual terms of the so-called Q matrix known to the matrix. 10 is also (Fisher's) information matrix and is defined by Q = A Σ 'A, where A is a theorem matrix and Σ is an error covariance matrix; vj is the actual term of v = ΑτΣΊηι and N defines x ^ n in the same neighborhood as (16) and ( The conditional distributions of other unknowns x2, X3, ... are similarly formed.

15

An important insight for the first aspect of the invention is that the conditional distribution of equation (18) is a piecewise Gaussian function of X [1]. (If Qn = 0, it is exponential.) The division into chunks comes from the last term, with the addition signs and the eigenvalues. Each of the eigenvalues has a break at xj = xi, because smaller values of xyn can be written to (xrXj), while larger values of x / .n the correct alternative notation is - (xj-xj) or (xj-xi). It is easy: to decompose the expression of equation (18) into separate parts and write it without any intrinsic value between each, where the positivity or negativity of (xi-xj) is known to all j; '; values in the neighborhood.

25] ··· _ Figure 2 illustrates schematically how a given arbitrarily chosen one-dimensional conditional distribution Dp for pixel • t of a two-dimensional rectangular lattice is piecewise Gaussian (or exponential). Pixel Close Neighborhood • ;;; is defined to include four other pixels, that is, those in contact with the · ·. ' 30 pixels per side. This yields four cut-off values, which are represented by vertical bars; as 201, 202, 203 and 204; at each border, the value of xf is equal to one of the other pixels. The boundaries divide the conditional distribution Dp into five parts 210,;, the parts of graph 210, the x-axis and the boundary line together delimiting the five two-dimensional zones 211, 212, 213, 214 and 215. Each of these belt zones has a well-defined area, at least assuming that the tails of the conditional and ejection positive and negative approach the x-axis sufficiently to obtain a finite area value corresponding to xj. n value range. According to statistical inversion theory, regularization utilizes probability distributions with finite surface area - one-dimensional conditional densities are the result of the use of such probability distributions, which also implies 5 finite integrals. In exceptional cases, it is also possible to set a positive and a negative limit on the "significant" values of Xj so that, on the other side of these limits, the conditional distribution receives only insignificant values; this ensures that finite area values are always obtained.

10 Creating a random value from a conditional distribution such as Figure 2 is a three-step process: 1) Calculate the area of each zone 211, 212, 213, 214, and 215. This happens quickly, since each band is part of a normal distribution (or an exponential 15 distribution), and the areas are readily obtained analytically or by a few simple calculations from the corresponding Gaussian (or exponential) cumulative probability tables.

2) One of the zones is selected at random but such that the probability of each zone 20 is proportional to its relative area. This can be done in many different ways. One very simple way is to associate each zone with a given region of real numbers between zero and one such that the extent of each region is proportional to the relative area of the corresponding zone and the regions: together fill the entire space, use a regular random number generator. : 25 to form a female number between zero and one and look at the area where the random number occurred.

3) Generating a random number from a selected range of possible xf values according to a probability density expressed by a normal (or exponential) function bounded to the corresponding range. This essentially corresponds to the known generation of random number ·· * cells from conditional distributions resulting from the use of ordinary Gaussian priorities (or exponential priorities), with the addition that this time the female number must be generated from a finite portion of the distribution. Depending on the shape of the zone, it may even be possible to ignore its Gaussian shape completely and approximate it, or more precisely, magnify it with a suitable, simpler function.

For example, it is possible to replace the small values of the normal distribution with a purely exponential distribution or to magnify a very even proportion of the normal distribution with an even distribution.

116524 12

The decomposition of a conditional distribution into piecewise Gaussian sections can thus combine the advantages of Gaussian priorities (or exponential or otherwise "easy" priorities) with the advantages of total variation type priorities. Priority guides computation toward solutions that are consistent with the structural features of many of the subjects being captured, which is characteristic of overall variation priorities. However, tedious numerical integrations of conditional distributions are unnecessary, whereby the reconstruction algorithm is much faster than corresponding prior art algorithms. Partitioning succeeds whenever the conditional distribution, when written in the form shown in equation (18), has absolute values of pixel values (or voxel values), absolute values of differences between pixel values (or voxel values) or higher values of pixel values (or voxel values).

Next, we analyze the theoretical part of the second aspect of the invention, which relates to the matrix computation needed in the image reconstruction algorithm. When the division Q1 of Q matrix 15 is computed in advance, one inner product and one addition of vectors must be calculated for each one-dimensional conditional distribution. One complete iteration thus requires a single matrix times vector multiplication with a relatively dense matrix: Q matrices are always much denser than their corresponding theoretical matrices (A matrices). The magnitude of this problem can be illustrated by an example where the tomography problem has 25,000 unknowns, 10,000 measurements, a theoretical matrix containing 1,106 nonzero elements and 80-106 nonzero elements. The Q matrix is symmetric, which means that only 40 · 106 nonzero elements are needed, but this also means; easily hundreds of megabytes. The number of processor operations is 2-80-106 FLOPS, or · 25,160 MFLOPS.

* * I

i '. * Another aspect of the invention is based on the realization that the A

* »I

; .. smaller than the Q matrix. This property of the theory matrix is a direct result of the sharp beam model essential to tomo-graph geometry. Figure 3 illustrates this fact in a two-dimensional form. Detector 301 consists of small spatially scattered detection elements. The quantum amount detected by one of the detection elements 302 is thought to be the result of a beam 303 passing along a straight line; Between the radiation source 304 and the detection element 302. On its journey, radius 303 passes through certain cells of the exaggerated coarse measuring grid 305 which are

::: 35 marked with a slash. Non - zero Effect on Theoretical Matrix

; : comes only from these shaded cells, which make up a small minority of all lattice points; hence the rarity of A. However, the Q matrix, which is essentially computed by ATA, is necessarily much more dense than A. The error covariance matrix joka that occurs in the complete definition of the Q matrix is diagonal in all practical cases and therefore has no effect on the rarity / density question.

The term actually needed to calculate one-dimensional conditional distributions, 5 is xTQx - (Ax) T (Ax). Instead of matrix Q, much less matrix A can be used.

Can be denoted A = (Ai; ...; A „), where the elements Ai; ... are a row of vectors. In this case, the calculation for the first coordinate is 10

Al {(A2 -, .. AAjxs ... xn)) = Ä [(Ax-Aixi) = Ä [(Ax) -Quxi. (19) This includes one vector times vector multiplication and basically one matrix times vector multiplication (Ax). In fact, the latter is not needed as the 15 new coordinates only need to update Ax. For example, in the case of the second coordinate,

Ax ^ Ax + A} (x ^ '- x ", d), (20) 20 or more generally

Ax -> A * + A,. (x ”K - xfd), ie (l, n), (21) • (» I *:; for the value of the current index i) Thus, one complete iteration requires two '*': 25 matrix times vector multiplication with matrix sparse A instead of dense <2.

In a practical algorithm, it may be necessary to recalculate the entire Ax from time to time to prevent the accumulation of rotation errors. A suitable frequency for recalculating Ax can be found by experimenting.

•: · 30 If you compare this with the example above, it is easy to see that the memory space '·; only 2106 nonzero items (non-zero A items and indices of those values) are required and the number of processor operations is only 4 ': MFLOPS.

; ; : 35 There is another way to speed up the computation in cases where the measuring grid is square ..: shaped and has an origin on the axis of rotation that defines the mutual motion of the tomographic imaging arrangement and the subject being photographed. In such a case, the theoretical matrix A calculated for the measuring direction a has substantially the same numerical content for the measuring directions α + π / 2, α + π and α + 3π / 2, taking into account only the relevant rotation. Therefore, the measurements need not be recorded for all the directions mentioned, but only for one direction. The calculation can be somewhat faster because four 5 times more calculations are performed for each element of A.

Figure 4 shows the main features of the tomographic imaging procedure. In step 401, the subject to be photographed is illuminated from multiple directions and corresponding measurement data is collected. In step 402, the reconstruction of the image begins by initializing the values of x to one of the initial values x0. In step 403, updated values of x are calculated, and in step 404 it is checked whether the end condition of the iterative process is satisfied. If not, return to step 403. The process loops through the loop of steps 403 and 404 until the end condition is met, after which the finished image is stored and / or displayed at step 405. The rendering step typically includes various sub-steps of rendering the described object two-dimensional objects from different angles. The steps shown in Figure 4 are general and common to all CT scans.

Fig. 5 illustrates how the first aspect of the present invention fits within the scope of Fig. 4. In a step of calculating updated x values, which in Fig. 4 is step 403, the steps of Fig. 5 are performed for each coordinate. In step 501, the formation of, ·, is piecemeal a Gaussian conditional distribution (or more generally, a conditional distribution which corresponds to some "easy" distribution in several discrete regions); ; to a point where limit values can be defined and the areas associated with each of them. In step 502, one area is randomly selected, taking into account the relative probabilities of the areas. In step 503:, · a random number is generated from the selected region, taking into account the probability density corresponding to a portion of the conditional distribution function in the selected region. The generated random number is stored in step 504 as the updated value of · 30 random variables that make up this coordinate.

Figure 6 shows the features according to another aspect of the present invention in the structure shown in Figure 5. In the step of generating a conditional distribution * function, which in Fig. 5 is step 501, the steps of Fig. 6 are performed. In step 601

35 is decided based on one of the pre-programmed decision criteria, whether it is time to recalculate input Ax or whether it is sufficient to update the previously calculated result. In the latter case, we proceed to step 602, whereby the previously calculated Ax is updated with respect to the last calculated new coordinate by multiplying the theoretical matrix A

116324 15 this line vector by subtracting the last calculated new coordinate value from its old value and adding the result to the previous Αχ value (see Equations 20 and 21). The updated Ajt value is then used according to step 603. An alternate route runs from step 601 to step 604 where a complete recalculation of the input Ax is performed, after which the result obtained is used as the new value of Ax in step 603.

Fig. 7 schematically shows an apparatus according to an embodiment of the invention. The device includes at least one controllable radiation source 701, at least one radiation detector 702, and a holder 703 for holding the object to be imaged relative to radiation source 701 and detector 702. The arrangement of the three parts shall be constructed in such a way as to allow for exposure from different angles. The basic principle can be implemented in practice using any of a variety of alternative methods. If there is only one source of radiation, either it or the target must be movable to irradiate the target from different directions; in this case, the device must have controllable moving means 704. The moving means may not be necessary if several radiation sources are located at different angles. The controllable radiation source 701, the radiation detector 702, the target holder 703 and the controllable moving means 704 may all represent prior art from prior art computed tomography systems.

20

Detector 702 is coupled to the input port of imaging computer 711. This data machine includes a measurement data acquisition interface 712, a data memory 713 for storing measurement data and calculated image data, and a program memory 714 for such computer programs; recording means for controlling an imaging computer 711 '/ 25 to reconstruct the image in accordance with the present invention. In addition, the imaging computer * 11 - * includes a database 715 for storing pre-computed data structures used in image reconstruction, such as theory matrices.

Of course, all of the aforementioned components of the imaging computer 711 are connected to at least one processor 716 that serves as the central control and processing of the computer · 30 units.

The operation of the imaging computer 711 is controlled through a user interface 717, which typically includes a keyboard 718 and a display 719. Between the processor 716 and the moving means 704, there must be a control interface 720 35 for control purposes for radiation source 701, radiation detector 702 and / or. :: to control the operation of the radiation source 701. The device may include various sensors 721 connected to a data acquisition interface 712 for monitoring environmental parameters.

116324 16

The difference with prior art computed tomography devices is, first, that the imaging computer 711 is programmed to perform the image reconstruction method of the present invention. However, the superior efficiency of the method of the invention over prior art methods 5 using statistical inversion also implies that the performance requirements set on the imaging computer 711 are much lower than hitherto considered mandatory. At the time of writing this disclosure, the present invention enables the use of a standard personal workstation to produce tomography images that are sufficiently accurate for clinical diagnostic use and yet require only a few seconds or few minutes to compute.

What the parts of the device of Figure 7 look like in practice depends largely on the subject being photographed. In medical imaging applications, the holder 703 is typically designed to assist the patient to remain motionless, while the controllable moving means 704 is arranged to move at least one of the radiation sources 701 or the detector 702. through a circular arrangement; alternatively, there are only one or at least a very small number of pairs of radiation source detectors, and the moving means are arranged to move the plug several times through the measuring arrangement and,, *. to rotate the log slightly each time about its longitudinal center axis. 'Ri. In industry or other applications that study the quality of workpieces, V! the holder 703 and the controllable moving means 704 may comprise complex robotic / 25 devices for automatically moving the radiation source 701, the detector 702, and the workpiece relative to one another. As a final example, satellite radiography may be used, in which the radiation sources 701 are placed in the satellites, or naturally occurring bodies in the outer space are used as the radiation sources. The detectors 702 are located in the satellites, whereby the controlled moving means 704 should be arranged to contain the satellite platform. In such a case, the generation and processing of the detection results involves a large amount of data transfer between the remote devices, and image reconstruction is clearly a separate process from data acquisition.

:; Figure 8 schematically shows a computer program product according to an embodiment of the invention. Process control means 801 is responsible for the smooth flow of the process stream by directing transitions between operating states, transmitting user-received commands to the appropriate operating units, monitoring the states of various subprocesses, implementing iteration criteria monitoring, etc. Exposure control means 812 is arranged to control possible movements 5 of radiation sources, detectors and target holders, as well as radiation from a controlled radiation source. Detection data acquisition means 813 is arranged to perform collection and storage of data obtained from radiation detectors. Theoretical matrix generating means 814 is arranged to assist the user in constructing the theoretical matrices. Priority forming means 815 is arranged to assist the user in setting priorities.

10

The conditional distribution generating means 816 is arranged to implement the principle of using a theoretical matrix instead of a Q matrix and updating the input Ax instead of recalculating it each time, as described in the second aspect of the present invention. The conditional distribution generating means 816 15 is further arranged to implement the principle that the conditional distribution is formed so as to consist of successive pieces which are "easy" to produce numbers. This part of the conditional distribution generating means 816 operates in accordance with the description of the first aspect of the present invention. Also, the area selection means 817 and the random number generating means 818 are arranged to operate in accordance with the description of the first aspect of the present invention, with proper selection of areas by "easy" conditional distribution. where and when generating random numbers from the areas so selected.

• ·. The image data storage means 819 is arranged to store the results of the finished image reconstruction structure and to prepare the image data so stored for display; ; ' under control of the interface means 811.

The examples given above should not be construed as limiting the applicability of the appended claims. For example, even the expression data; The collection process and the picture reconstruction process are described above as part of one, essentially uniform process, the invention does not exclude that detection data is collected in one location on one device and the collected data is sent to another location and / or another device for image reconstruction. The first and second aspects of the invention may also be distinguished; for example, if: 35 the object to be described is such that Gaussian priors can be used to obtain naturally normal conditional distributions, it is not necessary to form a conditional distribution consisting of successive fragments. However, in this case, another aspect of the invention, e.g. the principle of using the Q-116324 18 matrix instead of the theoretical matrix and updating the input Ax instead of recalculating it each time. On the other hand, if a new mathematical method is found to calculate the input Ax or to use the Q matrix efficiently, the application of another aspect of the invention may become unnecessary. In this case, the principle of avoiding difficult numerical integrations by the use of "easy" conditional distributions in a piecemeal fashion can significantly improve the overall efficiency of the method.

Another generalization of the previous descriptions concerns the mathematical method of finding the best possible solution for image data x. In the above description, 10 are focused on the so-called. conditional mean, or CM estimation, but other methods are known. The search for the best possible solution can be viewed as a very broad multidimensional optimization problem, which can be solved, for example, by using the so-called. MAP estimation (maximum a posteriori) accepting the x that maximizes the probability density function D (x \ m) over all 15 values of x. At the time of writing this description, CM estimation is preferred to MAP estimation because CM estimation has better error properties for the least squares method.

• * »» I t

Claims (24)

116324 A method for generating a tomographic image by statistical inversion, comprising the steps of: a) generating (403) a conditional distribution for each coordinate of an image to be formed, b) generating (403) a random value from the generated conditional distribution, and c) using (403, 405) a random value as a valid value for the pixel associated with the coordinate for which the conditional distribution was formed; 10 characterized in that - step a) comprises a sub-step of dividing (501; 201, 202, 203, 204) a conditional distribution (210) into parts, and - step b) comprising a sub-step of selecting (502) one of said parts, from which a random value is generated (503). 15 A method according to claim 1, characterized in that step a) comprises a sub-step of dividing (501; 201, 202, 203, 204) a conditional distribution into parts which themselves follow a normal distribution, whereby the conditional distribution (210) becomes piecemeal Gaussian. 20 A method according to claim 1, characterized in that step a) I 1 I *: includes a sub-step in which (501; 201, 202, 203, 204) the conditional distribution is divided into portions, each of which follows any of the following distributions: , normal distribution, exponential distribution. ; 25 A method according to claim 1, characterized in that step a) comprises a sub-step in which the outer portions include restrictions to ensure that the conditional distribution is limited. 5. The method of claim 4, further comprising scaling at least portions of the conditional distribution in the step of including constraints on the outermost portions to ensure that at its extreme values the conditional distribution approaches zero. Method according to claim 4, characterized in that, in the sub-step of including constraints on the outermost parts, the boundaries are determined for the minimum and maximum allowed value of the conditional distribution variable. 116324 A method according to claim 1, characterized in that step b) comprises sub-steps in which: - calculating the area of each zone (211, 212, 213, 214, 215) bounded by each part of the conditional distribution (210); sectors as part of the probability weighting of the components when selecting the part from which the index is generated. The method of claim 1, characterized in that step a) comprises the sub-steps of: 10. calculating a conditional distribution using a theoretical matrix that determines the proportion of measuring grid cells in the detection results obtained from the individual detection elements, and - at least in some cases said theory matrix jax is a vector of random variables whose value is to be determined in the image generation process by updating (602) said input Ax, which is performed by summing the previously obtained input Ax with a new calculated value of said line matrix A and said vector x and the product of the difference between the previously calculated value of said vector x. Method according to claim 8, characterized in that the method 20 comprises the step of recalculating (604) the entire product A x instead of merely updating when the touching condition (601) of such recalculation is fulfilled. »« Ψ • · * 10. A device for generating a tomography image, which device is arranged to utilize an I * »; ' ; a statistical inversion and the apparatus comprising: * * · 25 - means (714, 715, 716) for generating a conditional distribution for each coordinate of the image to be formed, - means (716) for generating a random value from the conditional distribution generated from the conditional distribution; for use as a valid value ·· 30 for the pixel associated with the coordinate for which the conditional distribution is "". formed; characterized in that the device comprises: - means (714, 715, 716) for dividing the conditional distribution into a plurality of parts, and ,,, · means (716) for selecting one of said parts as a part of which a random value. Apparatus according to claim 10, characterized in that said means (714, 715, 716) for dividing the conditional distribution into a plurality of portions is provided by a conditional distribution of the ground 116324 into portions which themselves follow a normal distribution, Device according to Claim 10, characterized in that said means 5 (714, 715, 716) for dividing the conditional distribution into a plurality of portions is arranged to divide the conditional distribution of the earth into portions, each of which itself follows any of the following distributions: Device according to Claim 10, characterized in that it comprises: 10. means (716) for calculating the area of each zone (211, 212, 213, 214, 215) bounded by each part of the conditional distribution (210) and means (716) to use the areas as probabilities weights of the parts when selecting the part from which the random number is generated. Apparatus according to claim 10, characterized in that: - means (714, 715, 716) for generating a conditional distribution are arranged to use a theoretical matrix which determines the proportion of measuring lattice cells to the detection results obtained from individual detection elements; for generating, there is provided a control-20 for replacing calculating an input Ax, where A is said theory matrix and x is a vector of random variables whose value must be determined in the imaging process by updating said input Ax, which is done by summing a 1, '; the aforesaid product Ax and the line vector from said theory matrix A and said *; ; know the new calculated value of vector x and the previously calculated value of said vector x! 25 is the product of the difference. Device according to Claim 14, characterized in that it is arranged: instead of merely updating, to recalculate the whole product Ax when the touching condition of such recalculation is fulfilled. 30 r: 1 | Device according to claim 10, characterized in that it is a medical imaging device comprising a controllable radiation source (701), a radiation detector; A '(702) and a holder (703) for holding the patient stationary during exposure and a controlled moving means (704) for controlling at least one of: a controllable radiation source (701), a radiation detector (702), and holder (703). Apparatus according to claim 10, characterized in that it is a tomographic log imaging device and includes an arrangement comprising a fixed radiation source (701) and a detector (702) 116324 and includes controllable moving means (703, 704) for moving the log axially to a fixed radiation source. with respect to the detector arrangement. A computer program computer program product for controlling a computer to generate a tomographic image by statistical inversion, said computer program product comprising: - program code means (816) arranged to direct the computer to form a conditional distribution for each coordinate of the image to be generated, 10. program code means (818) produce a random value from the generated conditional distribution, and program code means (819) configured to direct the computer to use the generated random value as a valid value for the pixel associated with the coordinate for which the conditional distribution was generated; Characterized in that the computer program product comprises: - program code means (816) arranged to direct the computer to divide the conditional distribution into a plurality of parts, and - program code means (817) arranged to direct the computer to select one of said parts as the portion from which the random value is generated. 20 A computer program product according to claim 18, characterized in that the program code means (816) arranged to control the computer forming a * * / '; * conditional distribution comprises program code means arranged to control::: the computer to divide the conditional distribution. to several parts which themselves:: 25 give a normal distribution, whereby the conditional distribution becomes piecemeal. · No. A computer program product according to claim 18, characterized in that the program code means (816) arranged to control the computer to form a conditional distribution comprises the program code means arranged to control the computer to divide the conditional distribution into a plurality of components, each - • 'The weather follows any of the following distributions: equal distribution, normal distribution, ex-. . * component distribution. *; 35 A computer program product according to claim 18, characterized in that the program code means (816) arranged to control the computer to form a conditional distribution of earth comprises program code means arranged to control the computer to include constraints on the outer portions to ensure a limited distribution. A computer program product according to claim 18, characterized in that the program code means (818) arranged to direct the computer to generate a random value from the generated conditional distribution comprises the program code means arranged to direct the computer to - calculate the area and sectors as probability weights of the components when selecting the part from which 10 random numbers are generated. A computer program product according to claim 18, characterized in that the program code means (816) arranged to direct the computer to form a conditional distribution comprises the program code means arranged to control the computer - to calculate the conditional distribution using a mathematical matrix that defines and - replacing, at least in some cases, computing the input Ax, where A is said theory matrix and x is a vector of random variables whose value must be determined in the image generation process, by updating said input Ax by summing the previously obtained input Ax and said is the row,. ·. the product of the difference between the new calculated value of the vector and said vector jc and the previously calculated value of said vector x. • * I A computer program product according to claim 23, characterized in that it comprises program code means (816) arranged instead of simply updating to recalculate the entire input Ax when the touching condition of such recalculation: .. is fulfilled. 116324