WO2008067842A1 - Method for generating at least two-dimensional tomographic image data, apparatus, computer program product - Google Patents

Abstract

Summarizing, the invention relates to a Method for generating at least two- dimensional tomographic image data of an object from a sinogram of the object, the method having the steps: - providing a sinogram s of an object, the object having an attenuation distribution map f(u, v); generating the attenuation distribution map f(u, v) of the object by minimizing the quadratic equation: F(f) = 1/2 fTAf - bTf using an iterative method, wherein A:= HTH and b:= HTs and wherein H is the projection matrix of the attenuation distribution map f(u, v) onto the sinogram s and HT is the backprojection matrix of the sinogram s onto the attenuation distribution map f(u, v) of the object; - deducing the at least two-dimensional tomographic image data from the attenuation distribution map, an apparatus and a computer program product.

Description

Description

Computer tomography allows the reconstruction of the three-dimensional inner structure of the investigated object from two-dimensional images showing different projections of the object. Applications range from medical imaging using x-rays to neutron tomography in various scientific and technical fields. Typically, tomographic image data are generated using a computer from a large amount of individual image data.

As far as computing is concerned, the problem is to provide sharp images rich in detail which are also easily discernible to the human eye. This demand is obvious in case of the clinic radiologist. In technical and scientific applications, such as materials research, for example, however, the output needs not only be qualitatively but also quantitatively correct.

Success of clinic treatment often depends on the rapid availability of reconstructions. The same demand often arises in industrial workflow and scientific procedure. Hence, computing time complexity and memory demand are also important figures of merit for any reconstruction algorithm. The memory constraint becomes ever more important because of the large amount of data produced by modern high resolution 3D imaging systems. Prior art

The most simple setup for tomography consists in a parallel beam penetrating an object which may be rotated about the vertical axis, assuming furthermore that the beam is attenuated by capture, not by scattering. In this case, horizontal "slices" of the sample do not influence one another and can be processed one after another. The following discussion is restricted to the reconstruction of one such slice.

The projection of the slice consists of a one-dimensional function. For any distance from the center of the beam it comprises the attenuation along that ray, which is, e.g., presented in figure 1 , showing two sinograms for different angles between the radiation and the object. In particular, figure 1 shows a cross-sectional view of a turbine blade as an exemplary object to be examined. The cross sectional view, as shown in figure 1 can also be a symbolic attenuation distribution map of the object in one plane, in particular as density profile of the object to be examined.

Put formally, for any angle φ there exists a function P, which maps the slice f(u, v) onto the projection P_φ(x). Here, u and v are the two-dimensional coordinates within the slice and x is the distance of the ray from the central beam through the axis or rotation.

Because attenuation can be considered an exponential process, using the logarithm of the projection turns the mapping into a linear function. This logarithm will be consequently referred to as the projection itself. The value of the projection at the point/distance x is then the integral of the attenuation coefficients along the ray with distance x from the central beam:

P_ψ(x) = / dz /( x cos ψ — z sin φ, x sin φ + z cos ψ )

= / άu J dv f(u, v) δ(x — u cos φ — v sin φ) This integral, known as the Radon transformation, could be inverted if it was known for any real value of φ and x. In practice, of course, only a finite number N of projections at different angles with a finite resolution of M pixels are recorded. The art of computer tomography consists in reconstructing f(u, v) as exactly as possible from this incomplete set of data.

As a result, f(u, v) will also be stored in a computer. The coordinates u and v are considered to be also of finite resolution U and V respectively. Exploiting the discreteness of all coordinates and the linearity of the projection mapping, the latter one may be written:

H : W U"V → _v mMN

With a sinogram s, i.e. the set projections for all different angles φ and the unknown attenuation coefficients f, the reconstruction problem can be stated in the form of a set of linear equations in matrix notation:

Hf = s

Naive methods for solving sets of linear equations cannot be used, however, because the system is large. Some high-resolution imaging equipment will provide projecting with a resolution of M ~ 10³. Assuming M - N - U - V, which is reasonable from a technical point of view and considering the Nyquist Condition, then H is a N² ® N² matrix. Gauss elimination would require O(N⁶) operations and O(N⁴) memory for storing the matrix, wherein H will most probably not be invertible.

In one approach, the projections are simply projected back onto the slice. The result is a very blurry image. This is the case, because any greyvalue collected by a ray travelling through a number of pixels of possibly different attenuation is now assigned to all of these pixels. Moreover, as the lines of projection are more dense in the center of the image, the image is artificially darker in the center, which is, e.g. shown in figure 2. In figure 2, there is shown a model of a turbine blade on the left hand side, as also shown in figure 1. On the right hand side of figure 2, there is shown the backprojection of the turbine blade.

In spite of its poor reconstruction properties, the merit of the backprojection is its speed and low memory demand. The following sniplets of pseudo code demonstrate how projection and backprojection can be carried out:

Projection (slice[N] [N] -» sinogram [N] [N]) V i, j: sinogram [i] [j] = 0 V j: co = Cos (π • j/N) si = Sin (π * j/N) V u:

V v: i = N/2+(u-N/2) . co+(v-N/2) • si sinogram [i] [j] += slice [u] [v]/N Backproj. (sinogram [N] [N] -» slice [N] [N]) V u, v: slice [u] [v] = 0

V j: co = Cos (π * j/N) si = Sin (π « j/N) V u: V v: i = N/2+(u-N/2) • co+(v-N/2) * si slice [u] [v] += sinogram [i] [j]/N

Obviously, projection and backprojection can be performed in O(N³) time with no more memory demand than the O(N²) space for the slice and the sinogram. For completeness and for a detailed understanding of the present invention, in the following, conventional tomographic methods are described.

Tomographic methods

A computer-tomography, also referred to as tomographic image data, typically is constructed by a two-step method. In a first step, measuring data are generated, which single elements do not comprise a direct three-dimensional information, taken alone. A single set of data does not code the information of a single point of a patient or an object. Rather, typically, it represents the information of points along a line through the patient/object. Only in a second step, a computer can reconstruct the entire three-dimensional structure of the patient/object, using said data. Depending on the specific process, the step of data generation can differ:

In medical tomographic devices of the first generation, as exemplary shown in figure 3, a radiation source 10, in particular an x-ray source 10 and perpendicular thereto a linear detector 12 rotate around an object/patient 14. The rays through the object/patient 14 are arranged in one plane, each. Each plane will be reconstructed, separately. After an entire rotation the object/patient 14 will be moved along the axis of rotation for a predetermined sequence. Therefore, slice by slice, an entire picture can be generated.

Using the positron-emission-tomography, a radioactive sample will be given to a patient, which sample decays by emitting positrons. Such a positron annihilates with an electron to two photons travelling under an angle of 180°. When localised on opposite sides of a detector 12, exemplary shown in figure 4, under coincidence, the positron source 10 is localised on the line connecting said coincidence detection point.

In parallel tomography, which includes neutron-tomography, the object will be rotated within a bundle of parallel rays, such that on the detector 12, as is also shown in figure 5, there are projected entire radiographies having different angles, each. A reconstruction will be generated in slices perpendicular to the axis of rotation.

The reconstruction of the tomography data conventionally is carried out by a computer. According to the method and the geometry of the data, the reconstruction methods can be different. However, in most cases, the problem can be simplified to parallel planes, which can be reconstructed independently and which are penetrated by rays.

In case of parallel tomography, these might be the planes perpendicular to the axis of rotation. Further to that, in case of medical tomography of the first generation, independently scanned planes are perpendicular to the axis of rotation. Hence, additionally, the rays have to be sorted according to the angle, since every record comprises a diversity of rays under different angles. There exist different methods, which can be used for reconstruction the entire three-dimensional volume within one step. Said methods are mentioned only briefly and shall not be further discussed. Rather, in view of the following, solely the cases with separate planes with parallel rays or a diversity of divergent rays shall be described.

Another possible way of grouping the possible reconstruction programs in order to describe a general method is based on the fact that the measured value of a ray can be described by an integral. Said integral can use, as a function of the integral, the quantity to be reconstructed, as it evolves along a ray. In case of the positron- emission-tomography each measured value is substantially proportional to an integral over the concentration of the radioactive marker. Mathematically, tomography is used to conclude a spatially dependent quantity μ from a line integral

Radon already examined this problem in 1917 and he also found a formula of inversion. However, practically applicable methods exist only, since the existence of powerful electronic data processing devices. Neutron Tomography

Neutron tomography is an example of parallel tomography. A neutron ray, which penetrates the object to be examined, is attenuated by absorption and scattering, as exemplarily shown in figure 6.

Along an infinite small distance ds along the ray, the intensity of the ray will be reduced by dl = - I μ ds, being dependent of the attenuation coefficient μ which has spatial dependency. After penetrating the sample, the ray has the remaining intensity

wherein I₀ is the intensity of the un-attenuated incoming ray. By measuring the intensity I₀ of the open ray without a sample, the line integral

can be obtained by normalizing and applying a logarithm. This is one of the line integrals, as mentioned above. By using a large variety of such line integrals and by tomographical methods, the attenuation coefficient μ can be obtained.

All ways of error avoidance and data correction, typically applied with regards to radiography, are applicable and useful for neutron tomography. This can be easily done, since the basis of tomography are complete radiographics and, since the necessary amount of line integrals is obtained by rotation of a sample within a ray around an axis of rotation. The axis of rotation is preferably perpendicular to the axis of the ray, which is schematically shown in figure 7. Due to the rotation of the sample, there are obtained conventional radiographics of the object under different angles of view, which can be seen in figure 7, upper right hand sight.

The thereby obtained data set is sorted, by summing the data referring to the respective horizontal detector lines. Each of the shown six data sets comprises the view of a horizontal layer of the sample under different angles (see figure 7, right hand side).

Each of these so called sinograms will be reconstructed tomographically, separately. Thereby, independent horizontal layers/slices are generated, which map the attenuation coefficient (see figure 7, lower right hand side).

By adding up the horizontal layers/slices to a three-dimensional rectangular parallelepiped, the therein contained information can be displayed by any possible views and within any possible slices, in a three-dimensional way (see figure 7, lower left hand side).

The Radon-Transformation

The coordinate systems for the two-dimensional problem have to be defined in a useful way. The fixed system will be fixed by two perpendicular axes in the horizontal plane. The respective coordinates are named u and v.

The description of the body reduces to a scalar function f of the real space, which does not necessarily have to be provided with an explicit physical meaning.

Regarding neutron tomography, it can be a spatially dependent absorption coefficient μ , as it is comprised in I - I₀ exp(- j//ώ) .

The function f does vanish identically outside the body, such that using a suitable axial scale, the carrier f in the space Ω is comprised in the interior of the unit circle S¹, as can be seen in figure 8. This example of the function f shows the so-called Shepp-Logan-Phantom, which is a standard test function for tomographic methods, conventionally known and used, which reminds on the interior of a skull.

The cross section of the object in the horizontal plane is tilted around the vertical y- axis, e.g. parallel to the paper plane. The bodily fixed u-v-system and the spatially fixed x-z-system are rotated against each other, which can be seen in figure 9. Further, line integrals of rays are regarded, wherein the line integrals are established by rotating the function f by the angle φ. Along a ray parallel to the z-axis having an intercept x on the x-axis, the line integral has the value

P_φ{χ) = $dzf(u(x,z,φ),v(x,z,φ)) . For each given angle φ there is obtained a function

P_φ having the argument x. Said function represents a projection of the function f under the angle φ onto the x-axis, which is schematically shown in figure 10.

Assuming such a projection for each angle φ, a new function P having the arguments φ and x is obtained using the function f having the arguments u and v. From the "ridge like" representation of f there is obtained an equivalent "ridge like" representation P, wherein the back of the ridge is obtained by the projections P_φ. For each single angle P (φ, x) = P_φ(x) which is schematically shown in figure 11.

The projection

R: f(u, v) → P(φ, x)

is the so called Radon-transformation. For simplicity, the Radon-transformed P(φ, x) is referred to as the projection of f(u, v), without confusing therewith the single projection under a specific angle. The Radon-transformed is also referred to as to sinogram and is typically displayed, as shown in figure 12. In such a representation as shown in figure 12, a point-like body transforms to a sinusoidal curve, therefore giving said representation said name. Since a sinogram is obtained by sinusoidal curves of different weight, amplitude and phase shift, using Fourier-methods, the original object could be recovered. An exemplary Fourier-reconstruction of a conventional Shepp-Logan-phantom is shown in figure 13.

The Inverse Problem

Radon showed that said transformation is reversible unique, such that from the projection P(φ, x) the original function f (u, v) can be reconstructed. However, said inverse problem is not trivial due to a plurality of reasons:

Firstly, the projection P(φ, x) is not actually known, i.e. not actually measured, as a function of the continuous variables φ and x, but will be sampled for a specific, discrete number of angles ψ_j, wherein also regarding x, a discretization on a limited number of pixels Xj is undertaken.

Secondly, the measured signal is necessarily erroneous, regarding neutron- tomography in particular due to noise according to the neutron statistics (or the lack of suitable neutron statistics), according to the basic noise of the detector and according to gamma-radiation.

Thirdly, deviations from the ideal geometry of the projection can occur. For instance due to scattered neutrons, distortion of the optics, imperfect alignment of the rotational axis of the object and the axis of the camera, etc. However, in case the real geometry is not coincident with the assumed geometry, the geometry of the objects cannot be reconstructed without error.

Finally, the typical tomographic reconstruction is comprised by the class of inverse problems, referred to as "slightly badly applied" problems. Badly applied problems have a high sensitivity regarding errors in the boundary conditions, the parameters of the model of inversion and the input data. In particular, said badly applied problems might be difficult to reconstruct, in case of limited angle/sub-angle tomography. In the following, there are represented direct reconstruction methods which are e.g. Fourier-transformation, the filterless backprojeciton and the filtered backprojection.

Direct Reconstruction

The Fourier-theory proves the invertibility of the Radon-transformation. Thereby, an explicit representation of the wanted function can be obtained, in which the backprojection meets the adjunct operation to the projection.

The backprojection can be used for direct reconstruction in real space. While the unfiltered backprojection, as a reconstruction method, provides visual and quantitatively undesirable results, the filtered backprojection, as a conventional method, avoids smearing of the image information by a suitable filtering of the sinogram before applying the backprojection.

The Fourier-Theory

The most direct reconstruction method follows from a proof of the invertibility of the Radon-transformation:

For the Radon-transformation, equivalently can be written:

P_φ(x) = J dz /( x cos φ — z sin φ, x sin φ + z cos φ ) = / du J dv f(u, υ) δ(x — u cos ψ — v sin φ)

Therefore, for any test-function g(x) the convolution characteristics apply: J dx P_φ(x) g(x) = J dx J du f dυ f(u, v) δ(x — u cos φ — v sin φ) g(x) — J du J dυ f(u, v) g(u cos φ + v sin ψ)

The convolution of relevance is the Fourier-transformed g(x) = (2π)^~U2 eχp(-ώr) , In this case, there holds:

JdX P_ψ{x) ^e^~ixr = JdU JdV /(«, V) -^e-^-^^vsin^^r

= V&ϊjdujdv /(u, υ) ^=_e-^iu(r∞B^ -fe._e-^iυl^rύnri

Obviously, the following Fourier-transformed are identical:

^v(^r) ⁼ V 2π/(r cos (/?, r sin ⁽^)

From the Fourier-transformed / a direct computation of / is possible.

Said provision however, cannot be applied for a practicable algorithm in a trivial way. The steps would be Fourier-transformation of P, calculation of / and, finally, the inverse Fourier-transformation of / to / .

The problem resides therein that / has to be given for the inverse transformation on cartesian grid points (in real life, x, z, r und φ are discrete). From P , the Fourier- transformed / is only known on a polar grid. The necessary interpolation provides small errors, which lead to visible artefacts in the inverse transformation, as can be seen in figure 13. The Concept of Backprojection

In case the direct application of Fourier-reconstruction is unavailable, it is possible to obtain useable approaches therefrom. Therefore, the wanted function f(u, v) will be explicitly written as

f{u, v) = (2π)-^χ Jda f db f{a, b) &"*^<*^>) a =: r cosfø) b =: r sinfø)

= (2Tr)-¹ Jdr frdφ f{rcos{φ), rsm(φ)) e^^ur∞s^^+υrsin^ o o

= (2π)-4 / dr / |r| d_¥> (2π)-¹/2 _jp_v>(_r) _eir(_ϋcos(_v»)_+υsin(_v))

0

= /dv_? (2π)"¹/² / dr FP_φ(r) _e*K^»«^«(v⁾+^»Bin( _v ⁾⁾ o

= J dφ F P_φ(x) x := u cos(<£>) + v sva(φ) o

= P^τ F P_φ(x)

With the filter F, which has the following properties in frequency space

and the backprojection

2π

P^T9{φ, r) = / dψ g(φ, u cos(9?) + v sάn(φ))

The backprojection transfers the sinogram having the coordinates (x, φ) into the bodily fixed system (u, v). As already shown with regard of the projection, the effect can be visualized for single, fixed angles φ, as schematically shown in figure 14.

For each fixed x the line integral integrated by the projection along u cos(φ) + v sin(φ) = x is reflected on said straight line by backprojection. The entire backprojection thereby comprises integration of all such backprojections in part for all angles φ (or cp_j, respectively).

Unfiltered Backprojection

The easiest of all reconstruction methods is the backprojection of the unfiltered sinograms. In the assumption

f(u, v) = P^T F P_ψ(x) w P^τ P_ψ{x) = / dip P_ψ{x = u cosfa) + υ sin (φ))

0

the wanted function f can be assumed directly from the measured data Pφ, without the detour over the frequency space.

As can be seen later, the measured data are only available for a limited number of angles ψ_j with j e {1 , ... N_j} and only for a limited number of detector elements Xj with i e {1. ... N,}.

Accordingly, the result f(u, v) will be stored as a matrix having a limited number of coordinates u_m and V_n with m e {1 , ... N_m} and n e {1 , ... N_n). Instead of the mapping x = u cos (φ) + v sin (φ) there has to be provided a respective weighting of the overlapping small areas.

Regardless of that, the discretized backprojection can be written

The stepwise summation of projections provides a successive and approximate reconstruction, as can be seen in figure 15, which shows reconstructions having 1 , 2, 4, 8, 32 and 128 projections.

The Filtered Backproiection

The simple assumption of the unfiltered backprojection, f(u, v) = P^τFP_φ(x) « P^τP_φ(x), already shows the object to be reconstructed. The reconstruction is blurry and, compared to the rim, too dark in the interior.

Filtering of the sinograms before the backprojection eliminates said defects. The necessary filter properties in frequency space

Fg(<P, r) = ttg(φ, r)

is provided by the Riesz-potential 7¹ and F = (^)^"1/¹. However, this amplifies the errors in the upper frequency range of the data, which are particularly affected by the discretization, as can be seen below. Accordingly, the data have to be additionally filtered. Actually, said filtering is supplied as convolution in real space. The necessary core

depends on the choice of the low pass F (r). For the ideal low pass filter F = Θ(π -\r\) the Ram-Lak-filter is obtained in discrete form:

according to Ramachandran and Lakshminarayanan. Fewer artefacts are produced by the filter

obtained by Shepp and Logan from F =

-smc(r/2) . The actual filtering occurs according to the convolution

(FPUx) = ∑ F(Ax) P_φ(x + Ax)

Ax

The result of the backprojection of the filtered sinograms is schematically shown in figure 16, which shows reconstructions using 1 , 2, 4, 8, 32 and 128 Shepp-Logan- filtered projections.

Typically a discretization of a projection can be performed to obtain a linear system.

--L <" I 5-SΣ Ni - 2 N_j

Object of the Invention

It is an object of the invention, to provide at least two-dimensional tomographic image data more efficiently, i.e. with the need for fewer computer resources and being more flexible with regards to the peripheral boundary conditions, such as room available, radiation available, etc. Said object is solved by the method according to claim 1 , the apparatus according to claim 15, the use of an apparatus according to claim 17 and the computer program product according to claim 18. Preferred embodiments are subject of the dependent sub claims.

Particularly, a discretization can be carried out, as described below.

Method according to an Aspect of the Invention

One aspect of the invention relates to a method for generating at least two- dimensional tomographic image data of an object from a sinogram of the object, the method having the steps:

providing a sinogram s of an object, the object having an attenuation distribution map f(u, v);

generating the attenuation distribution map f(u, v) of the object by minimizing the quadratic equation:

F(f) = 1/2f^TAf - b^Tf

using an iterative method, wherein

A: = H^TH

and

b: = H^τs

and wherein

H is the projection matrix of the attenuation distribution map f(u, v) onto the sinogram s and

H^τ is the backprojection matrix of the sinogram s onto the attenuation distribution map f(u, v) of the object;

deducing the at least two-dimensional tomographic image data from the attenuation distribution map.

A sinogram s comprises one or more intensity data sets. Each intensity data set can be, e.g., a vector which comprises one or more intensity values. The intensity values can be intensity values, which are generated, when radiation is penetrating the object. Each intensity data set corresponds to a specific radiation angle. In other words, radiation penetrating the object will be attenuated due to the attenuation of the object, which is described by the attenuation distribution map f(u, v). After penetrating the object, the attenuated radiation, which is attenuated according to the attenuation distribution map f(u, v) of the object, is detected by a detector. According to the attenuation of the radiation, intensity of the attenuated radiation detected by the detector can vary. Given an origin of a coordinate system, which, e.g. corresponds to the origin of the detector, the intensity data set represents the intensity of the attenuated radiation detected by the detector with respect to the origin.

Also, each intensity data set is generated under a specific radiation direction, wherein the radiation direction represents the relative angle between the object and the ray along which the radiation travels. For example, the alignment of the rays with respect to the detector can be fixed and the object can be rotated with respect to the incoming radiation, thereby changing the radiation angle, i.e. the radiation direction. In other words, at first instance, a coordinate system of the object and a coordinate system of the system (comprising the detector and the incoming radiation, which has a specific direction with respect to the detector) can be aligned. In a next step, the object can be rotated such that the coordinate system of the object and the coordinate system of the system are arranged having a specific angle with respect to each other. A further intensity data set can be established from said rotated object. In the coordinate system of the object, the intensity distribution map f(u, v) is constant for any rotational angle of the object (see e.g. figure 1 , figure 8, figure 9, figure 10 and figure 11 ). Said rotation substantially represents a change of a relative angle between the object and the system comprising radiation source and detector. Accordingly, said rotation also changes the direction of the radiation with respect to the object.

Accordingly, each intensity data set shows the attenuation of radiation penetrating the object in one plane (under a specific radiation direction). The radiation will be attenuated according to the attenuation distribution map f(u, v). The attenuation distribution map f(u, v) can be considered a cross-sectional view of the attenuation of the object. In the simplest case, the attenuation distribution map f(u, v) can be a density distribution map. However, further parameters can influence the attenuation distribution map f(u, v) and can be incorporated in the attenuation distribution map f(u, v), such as scattering of the radiation, etc. In other words, the attenuation distribution map f(u, v) for the same object can change depending on the radiation penetrating the object. For example, the attenuation distribution map f(u, v) of the object regarding x-ray radiation for generating the intensity data set can be different from the attenuation distribution map f(u, v) of the object, when generating the intensity distribution set using neutron radiation.

Therefore, the sinogram s preferably comprises two or more intensity distribution sets, wherein the radiation directions (which can be another expression for the alignment of the coordinate system of the object and the coordinate system of the system) differ from each other for each intensity data set.

Hence, preferably, each intensity data set is a one-dimensional data set, comprising the intensity of the radiation after penetrating the object as a function of the distance from an origin of the detector. The attenuation distribution map f(u, v) represents a two-dimensional data set comprising the attenuation factor of the object as a function of the coordinates of the object in the coordinate system of the object. A respective coordinate system can, e.g., be a cartesian coordinate system, a polar coordinate system, etc.

Furthermore, the step of providing a sinogram s can comprise shining radiation onto an object and detecting the radiation after penetrating the object, i.e. detecting the intensity of the radiation after penetrating the object. Alternatively/additionally, a sinogram data can be at least partly taken from a data base. In particular, there can be provided a specific step of generating a sinogram and a further step of generating the attenuation distribution map from the sinogram using the quadratic equation:

F(f) = 1/2f^TAf - b^Tf.

Moreover, the at least two-dimensional tomographic image data, in the simplest case, can be identical to the attenuation distribution map f(u, v). However, the two- dimensional tomographic image data can differ from the attenuation distribution map f(u, v), e.g. by introducing false colours, weighting of the displayed data, etc.

According to the present invention, it was found out, that, regarding in particular filtered backprojection, an alternative approach for deducing the at least two- dimensional tomographic image data can be applied. According to the invention, it was in particular found out that the naive and simple backprojection can be turned into an accurate reconstruction algorithm, advantageously still sharing the time and memory complexity of the backprojection. In particular, it was found out that starting from the simple description of the projection:

Hf = S

with H: R^uv → R^MN applying the backprojection H^τ to both sides of the equation, the composite operator H^TH is a square metric (it maps R^uv onto itself), symmetric by construction (H^TH)^T = H^TH) and positive-definite (f^TH^THf > 0).

Further, according to the invention a quadratic term F(f) = 1/2f^TAf - b^τf has a gradient F¹ (f) = 1/2 A^τf + 1/2 Af - b which vanishes at the minimum of F. Considering A being symmetric, there holds 0 = 1/2 A^τf + 1/2 Af-b = Af-b. Therefore, minimizing the quadratic form is equivalent to solving the set of linear equations of F(f).

Further, according to the invention, it was found out that the above indicated equation (H^THf = H^τs) can be simplified by defining:

A: = H^TH

and

b: = H 's.

Therefore, according to the invention it was found out, instead of deducing the attenuation distribution map f(u, v) by applying the backprojection and/or filtered back projection, etc., a different approach is realized. In particular, instead of solving backprojection and/or filtered backprojection directly, a supplementary equation, namely F(f) = 1/2 f^TAf - b^τf is minimized.

In other words, according to the present invention, it was found out that the attenuation distribution map f(u, v) can be obtained not only directly, as conventionally known. Rather, it was found out that a supplementary equation can be established and by minimizing said supplementary equation, the attenuation distribution map f(u, v) and/or the tomographic image data can be obtained.

In order to minimize the above supplementary equation, i.e. the above displayed quadratic form, a variety of iterative methods can possibly be applied.

According to a preferred embodiment of the present invention, the iterative method uses conjugate gradients and/or conjugate directions and/or a steepest decent, etc.

In other words, methods to minimize the above indicated quadratic form are, e.g., the method of steepest decent, the method of conjugate directions, the method of conjugate gradients, and so forth. However, it was found out that, advantageously, the method of conjugate gradients shows the best performance.

Using pseudo code, the preferred method for generating at least two-dimensional tomographic image data can be described in the following way:

i = 0 (initialisation) f = b (or f = 0 or ...) r = b - Af (the residuum d = r (initial direction)

Onew ⁼ r r i = 1 , 2, 3 ... (iteration) q = Ad α = δ_new/(d q) (optimal step size) f = f + αd (step to new point) r = r - αq (the new residuum)

Oold ⁼ O_new s- _ T Onew — I I β = δnew/δold d = r + βd (the new direction)

In other words, the input of the conjugated gradient algorithm is the unfiltered backprojection b which can be advantageously computed in time O(N³) and memory O(N²).

The relevant operation in the algorithm is the matrix product Ad. The matrix product Ad can be established by simply applying the projection and the backprojection on d. Advantageously, neither the sparse matrices H and H^τ nor the product H^TH need actually be computed or stored. Further advantageously, the remaining scalar products and vector additions can all be performed in a small or neglectable amount of time, in particular in an amount of time O(N²).

The number N in that case, relates to the number of intensity data sets each sinogram s comprises. In other words, N relates to the number of intensity data sets obtained under different radiation angle/rotational angle of the object relative to the radiation direction. According to the above and following terminology, N = N_j.

Moreover, advantageously, the complexity in computing time is thus T_CG ~ (3 + 2C)T_Bp, where T_Bp is the computing time for the naive backprojection and C is the number of iterations. The latter one depends only on the desired image quality and does not influence asymptotic behaviour. In practice, C is of the order of 10.

According to a further preferred embodiment, the sinogram s comprises at least two intensity data sets, wherein each intensity data set is generated from radiation penetrating the object under specific radiation direction with respect to the object and wherein for all intensity data sets generated from radiation having respective radiation directions φ_k, φι with k ≠ I, there is provided:

φ_k- φι < 180°, preferably

φ_k- φι < 120°, further preferably

φ_k- φι < 90°, most preferably

Φk- φi < 10°.

For example, all the intensity data sets can be obtained with an angular range within one interval, e.g. the angular range can be a continuous set of (discrete) angles between approximately 20° and approximately 120°. In particular, it is possible that the intensity data sets are not generated from radiation lying within one direction/angle interval. Rather, it can be possible that a first part of the intensity data sets is generated from radiation within a first interval and a second number of intensity data sets is generated from radiation within a second direction/angle interval which can be different from the first interval.

According to conventional tomographic methods, the object typically has to be rotated by 360° respect to a fixed radiation source, i.e. within the radiation, (and/or, as indicated above, the radiation source and the detector can be rotated with respect to the object to be examined). According to a preferred embodiment of the present invention, advantageously, it is possible to provide a sinogram by rotation with an angle of rotation not greater/less than 180°, preferably not greater/less than approximately 120°, further preferably not greater/less than approximately 90°, most preferably not greater/less than approximately 10°.

In other words, the coordinate system of the object can be arranged with respect to a fixed coordinate system of the system comprising radiation source and detector under a specific angle. Radiation is applied to the object and an intensity data set can be generated. Next, the object can be rotated, wherein the difference between the initial angle between a specific axis of the coordinate system of the object and the respective axis of the coordinate system of the system (comprising radiation source and detector) is changed. Said new angle between said two axis can be referred to as the radiation angle or the radiation direction. With respect to figures 8 and 9, in figure 8, the coordinate system of the body is identical to the coordinate system of the system. In figure 9, the body is rotated, the angle between a v-axis of the coordinate system of the body and the z-axis of the coordinate system of the system is referred to as φ, i.e. can be referred to as radiation direction, radiation angle, etc. Therefore, in figure 8 the radiation angle/direction can be defined a 0°.

On the other hand, with respect to figure 1 , the radiation angle can be the angle of an axis of the coordinate system of the body (not shown), such as the v-axis (as shown in figure 8 and 9) and the rays of incoming radiation, represented by the arrows, in figure 1. By rotation of the object, as shown on the right hand side in figure 1 , the angle φ between the v-axis (not shown) and the arrows representing the radiation, is changed. Therefore, also the radiation direction (with respect to the object) changes. The same holds, when rotating/moving the radiation source and/or the detector.

In particular, when changing the radiation direction, the object and/or the set comprising the radiation source and the detector are preferably rotated around an axis which is perpendicular to the plane comprising the u-axis and the v-axis of the attenuation distribution map.

Moreover, according to a preferred embodiment of the present invention, instead of providing intensity data sets for an entire/complete angular range between, e.g. 0° and 360°, it is possible to provide intensity data sets within two or more intervals, e.g. to provide intensity data sets from the interval of angular directions between 20° to 40° and to provide intensity data sets within the interval of 90° and 120°. This might be the case, since between e.g. the interval of 0° and 10° the object cannot be penetrated by radiation and/or the penetration of the object by radiation is bad. Accordingly, within the interval between 40° and 90°, the radiation penetrating the object with radiation might not be possible and/or bad. In particular, this could be due to mechanical reasons, such as that the object and/or the radiation source cannot be rotated entirely around the object. On the other hand, this might have e.g. medical reasons, such that specific parts of the body should not be exposed to radiation.

However, according to the present invention, even under a small interval of radiation directions, advantageously, the attenuation distribution map f(u, v) can be obtained sufficiently well. Moreover, even without a continuous interval of radiation directions, in particular when providing a sinogram having intensity data sets from discontinuous radiation direction intervals, the attenuation distribution map f(u, v) can be obtained, sufficiently.

According to a further preferred embodiment, for calculating the iterative method, specific trigonometric calculation operations are used, which cancel the typical rounding errors of operating microprocessors.

Particularly preferably, a wrap around operation is carried out on the cos(x) and the sin(x) functions of the standard libraries.

Examplarly, e.g. for angles φ^ π*— , such a wrap around operation can be

N described in pseudo code:

if (6*j == 0*N) co = 1.0 else if (6^*j == 1*N) co = 0.866025403784438646763723170753 else if (6^*j == 2^*N) co = 0.5 else if (6^*j == 3*N) co = 0.0 else if (6^*j == 4^*N) co = -0.5 else if (6^*j == 5*N) co = -0.866025403784438646763723170753

else co = Cos(π*j/N)

According to a further preferred embodiment, the resolution of the attenuation distribution map f(u, v) is virtually enhanced when projecting the attenuation distribution map f(u, v) onto the sinogram s.

Most preferably, each point of the attenuation distribution map f(u, v) is equally divided into a subset of virtual points and each virtual point is projected onto the sinogram s, wherein the weights for the virtual points are equal.

In other words, each of the elements of the slice of the object which will be projected is covered by an equally distributed point grid. Each point of the equally distributed point grid is projected having a specific weight. Each weight has the value of one divided by the number of points. Since said weight is applied onto the entire projection, to each detector element corresponds a natural number of points. According to a preferred embodiment of the invention, the projection matrix H and/or the backprojection matrix H^τ are compressed projection matrices.

Most preferably, the projection matrix H comprises entries solely for elements of the sinogram s, onto which points of the attenuation distribution map are projected.

Accordingly, also the backprojection matrix H^τ comprises individual values, only.

In particular, preferably, when calculating a projection (backprojection accordingly) within one loop all points of a picture/image are projected onto an intensity data set, for one radiation direction/radiation angle. Accordingly, within one loop, all points of a picture/image are projected on the sinogram s, for all direction angles. The coordinates in the sinogram represent the one, the coordinates in the picture the second index of the projection matrix. Therefore, exactly the non-furnishing entries of the matrix are specified. In an upfront calculation these coordinates can be stored in the following way:

loop over all angles φ - loop over the image points, i.e. the points of the attenuation density map f(u, v) project each point of the attenuation distribution map f(u, v) under the angle φ onto (plural) entry (entries) i in the sinogram; store the smallest i; - store for each consecutive i the number of projected points; store a zero.

In particular, the above method can apply, when enhancing the resolution, as described above, i.e. by covering each point to be projected by a grid of a predeterminable number of grid points. For example, such a number can be the square of any prime number (i.e. a square having a length of a side equal to a prime number). In particular such a number can be 29 x 29 grid points. In other words, when artificially enhancing the resolution, each surface element, which has to be projected, will be represented, e.g., by 29 x 29 points, 31 x 31 points, 37 x 37 points etc., which are projected individually. Hence, preferably, each picture/image element (i.e. each element of the attenuation density map f(u, v)) will be mapped onto the sinogram weighted according to its area and vice versa.

Preferably, a point radiation source or a radiation source radiating parallel radiation is used.

Particularly preferably, the radiation is neutron radiation or positron radiation or electron radiation or proton radiation or ultrasound or a combination thereof.

In other words, preferably, any known and/or conventionally used radiation can be used. In particular, conventionally known positron emission tomography can be applied using a method according to a preferred embodiment of the invention.

More preferably, the radiation is electromagnetic radiation having a wavelength from approximately 5 femtometres to approximately 10 nanometres. In other words, the electromagnetic radiation can be x-ray radiation, particularly soft x-rays or hard x- rays.

Preferably, the direction of the radiation relative to the object is changed by rotating the object and/or rotating the radiation source.

Preferably, the radiation source can be rotated around a center of rotation which is outside of the rotation source. Additionally/alternatively, the radiation source can be rotated around a center of rotation, which lies within the radiation source.

According to a further preferred embodiment, the tomographic image data are three- dimensional data, in particular generated from a plurality of attenuation distribution maps f(u, v). Apparatus according to an aspect of the invention

Another aspect of the invention relates to an apparatus for generating at least two- dimensional tomographic image data of an object from a sinogram s of the object, comprising:

a sinogram provision means adapted for providing a sinogram s of an object, the object having an attenuation distribution map f(u, v);

an attenuation distribution map generating means adapted for generating the attenuation distribution map f(u, v) of the object by minimizing the quadratic equation:

F(f) = V₂ f^TAf - b^τf

using an iterative method, wherein

A:= H^TH

and

b:= H^τs

and wherein

H is the projection matrix of the attenuation distribution map f(u, v) onto the sinogram s and

H^τ is the backprojection matrix of the sinogram s onto the attenuation distribution map f(u, v) of the object and a tomographic image data deducing means for deducing the tomographic image data from the attenuation distribution map.

As preferred embodiment(s) of the invention, the specific elements of the apparatus can comprise conventional elements known in conventional computer technology and/or conventional tomography technology. In particular, the attenuation distribution map generating means can comprise a microprocessor and/or respective controllers and bus systems. Moreover, the sinogram provision means can comprise a conventional radiation detector, such as a CCD-chip, etc. and/or a conventional optical drive/reader and/or optical storage medium, a conventional magnetic drive/reader and/or magentic storage medium, a conventional network drive and/or network connection, etc. Furthermore, the tomographic image data deducing means can comprise a microprocessor and/or respective controllers and bus systems.

Preferably, the sinogram provision means comprises a radiation source and/or a radiation detector.

Preferably, the attenuation distribution map generating means is adapted to use conjugate gradients and/or conjugate directions and/or steepest decent calculation method(s).

Preferably, the sinogram provision means is adapted for providing a sinogram s comprising at least two intensity data sets, wherein each intensity data set is generated from radiation penetrating the object under a specific radiation direction with respect to the object and wherein for all intensity data sets generated from radiation have been respective radiation directions ς>_k, φι with k ≠ l, there is provided:

(p_k- φι < 180°, preferably

φ_k- φι < 120°, most preferably φk - φι < 10°.

Further preferably, the attenuation distribution map generating means is adapted for calculating the iterative method using specific trigonometric calculation operations which cancel the typical rounding errors of operating microprocessors.

Most preferably, the attenuation distribution map generating means is adapted for carrying out a wrap around operation on the cos(x) and the sin(x) functions of the standard libraries.

More preferably, the attenuation distribution map generating means is adapted for virtually enhancing the resolution of the attenuation distribution map f(u, v), when projecting the attenuation distribution map f(u, v) onto the sinogram s.

Particularly preferably, the attenuation distribution map generating means is adapted for equally dividing each point of the attenuation distribution map f(u, v) into a subset of virtual points wherein each virtual point is projected onto the sinogram s and wherein the weights for all the virtual points are equal.

Further preferably, the attenuation distribution map generating means is adapted for compressing the projection matrix H and/or the backprojection matrix H^τ such that said matrices are compressed projection matrices.

Particularly preferably, the projection matrix H comprises entries solely for elements of the sinogram s, onto which points of the density distribution map are projected.

According to a further preferred embodiment, the apparatus comprises a point radiation source or a radiation source radiating parallel radiation.

More preferably, the radiation source is a neutron radiation source and/or a positron radiation source and/or a electron radiation source and/or a proton radiation source and/or an ultrasound radiation source. Most preferably, the apparatus can be adapted to provide conventionally known positron emission tomography.

Further preferably, the radiation source is an electromagnetic radiation source, having a wavelength from approximately 5 femtometres to approximately 10 nanometres.

According to a further preferred embodiment, the apparatus is adapted such that the direction of the radiation relative to the object is changed by rotating the object and/or by rotating the radiation source.

Further preferably, the tomographic image data deducing means is adapted to deduce three-dimensional tomographic image data in particular from the one or more attenuation distribution maps f(u, v).

Use of an apparatus according to an aspect of the invention

A further aspect of the present invention relates to the use of an apparatus according to the present invention in a security system, in particular in an airport security system.

Advantageously, in airports, at customs, in trains stations, etc. an easy and reliable way of monitoring objects, such as luggage, cameras, computers, etc. can be provided. The luggage can be checked very quickly, since only a fraction of a rotation angle, e.g. such as 120°, 90°, 60°, 45°, 30°, 10° can be used to obtain a complete image of the object monitored. Therefore, advantageously, the object can monitored in a fraction of the conventionally necessary time. Moreover, three-dimensional image data can produced instead of the conventionally possible two-dimensional image data. This is in particular the case, when a conveyor belt carries items through an x-ray machine, preferably without the need for an actual rotation of the object. Preferably, intensity data is acquired in fan-beam-geometry exploiting the translatory motion of the conveyor belt in front of the x-ray source. Advantageously, tomographic image data, in particular, three-dimensional tomographic image data can be reconstructed from incomplete input data. Therefore, the apparatus according to the present invention can be an important tool to deter risks in aviation, in particular in the detection of explosives in carry-on items going through the x-ray system or the detection of explosives in air cargo.

Moreover, e.g. at customs, advantageously, large objects, such as cars, trucks, etc. can be monitored in a fast and efficient way, allowing to produce three-dimensional image data. Therefore

Computer program product according to an aspect of the invention

A further aspect of the invention relates to a computer program product, which, when loaded in the memory of a computer and executed by the computer, performs the method according to the invention.

Preferred embodiments of the present invention are exemplarily described regarding the following figures. Although specific embodiments can be described independently, individual features and functionalities of the respective embodiments can be combined without prejudice, in order to provide further preferred embodiments.

Fig. 1 : shows a schematic view of an examination of an object under different radiation directions/angles;

Fig. 2: shows a schematic view of an object and a schematic view of a backprojection;

Fig. 3: shows a schematic view of an apparatus; Fig. 4: shows a schematic view of an apparatus ;

Fig. 5: shows a schematic view of an apparatus;

Fig. 6: shows a schematic view of a Shepp-Logan-Phantom and the attenuation of radiation;

Fig. 7: shows a schematic view of neutron tomography;

Fig. 8: shows a schematic view of a Shepp-Logan-Phantom and respective coordinate systems;

Fig. 9: shows a schematic view according to figure 8, wherein the Shepp- Logan-Phantom is rotated;

Fig. 10: shows a schematic view of a Shepp-Logan-Phantom and a corresponding intensity data set;

Fig. 11 : shows a plurality of Shepp-Logan-Phantoms and corresponding intensity distribution sets;

Fig. 12: shows a schematic view of a Radon-transformation;

Fig. 13: shows a schematic view of a direct Fourier-reconstruction;

Fig. 14: shows a schematic view of simple backprojection;

Fig. 15: shows a schematic view of unfiltered backprojection having several projections under different angles;

Fig. 16: shows a schematic view of filtered backprojection from a plurality of Shepp-Logan filter projections;

Fig. 17: shows a schematic, exemplary view of a possible discretization of the coordinate system;

Fig. 18: shows a schematic, exemplary view of a matrix;

Fig. 19: shows a schematic, exemplary view of a matrix product;

Fig. 2OA: shows a schematic, exemplary view of an iteration using a conjugate gradient method;

Fig. 2OB: shows a schematic, exemplary view of another iteration using a conjugate gradient method;

Fig. 2OC: shows a schematic, exemplary view of another iteration using a conjugate gradient method;

Fig. 21 : shows a schematic, exemplary view of a difference of two reconstructed tomographic image data;

Fig. 22: shows a schematic, exemplary view of an error as a function of the interation step;

Fig. 23: . shows a schematic, exemplary view of tomographic image data using filtered backprojection;

Fig. 24: shows a schematic, exemplary view of tomographic image data using a method according to a preferred variant of the present invention;

Fig. 25: shows a schematic, exemplary view of an error as a function of the interation step;

Fig. 26A: shows a schematic, exemplary view of tomographic image data using a method according to a preferred variant of the present invention or using unfiltered backprojection;

Fig. 26B: shows a schematic, exemplary view of tomographic image data using a method according to a preferred variant of the present invention or using unfiltered backprojection;

Fig. 26C: shows a schematic, exemplary view of tomographic image data using a method according to a preferred variant of the present invention or using unfiltered backprojection;

Fig. 27A: shows a schematic, exemplary view of devices according to preferred embodiments of the present invention;

Fig. 27B: shows a schematic, exemplary view of a device according to preferred embodiment of the present invention at different time stages;

Fig. 28: shows a schematic, exemplary view of a way to discretize fan beam geometry;

Fig. 29: shows a schematic, exemplary view of tomographic image data;

Fig. 30: shows a schematic, exemplary view of tomographic image data;

Fig. 31 : shows a schematic, exemplary view of weighted projection;

Fig. 32: shows a schematic, exemplary view of a reconstruction of noisy data; Fig. 33 shows a schematic, exemplary view of tomographic image data;

Figs. 34A to 34C: show a flow diagram of a method according to a preferred variant of the invention;

Fig. 35: an exemplary system comprised in an apparatus according to a preferred embodiment of the invention.

Detailed description of the Figures

The discretization of the projection which is exemplary shown in figure 17 (displaying a line of a detector 12) describes projection and backprojection as linear functions, which can be described as matrix operations on vectors. In particular, figure 17 shows a detector 12 being an exemplary line detector. Further to that, a grid 16 is displayed. Additionally, a projection is symbolised by an arrow 18. The grid 16 can represent the attenuation distribution map f(u, v), wherein each entry of the grid can represent an element of the attenuation distribution map f(u, v). Hence, the grey shaded areas, connected by the arrow 18 can represent a projection of the member of the attenuation distribution map onto the line detector 12. Further, there is displayed the coordinate system of the attenuation distribution map, made of axis u and v and also the radiation angle/direction φ of the radiation to be detected. The radiation angle/direction φ can also represent the plane of the attenuation distribution map, i.e. the cross section of the body, with respect to the a coordinate system of the detector and/or the radiation source.

Following that, according to a preferred embodiment of the invention, tomography can be formulated as a linear problem. Preferably formulated as quadratic, positive form, there is provided the algebraic method according to a preferred embodiment of the present invention, which establishes the reconstruction iteratively using e.g. the method of conjugate gradients. In contrast to filtered backprojection, advantageously this method is suitable, when only incomplete information is provided, such as measured data, which is provided for one, two, three, etc. (i.e. a limited number of) angular intervals only. Further advantageously, said method is suitable for a fan beam geometry of non-continuous angular intervals. Regarding incomplete information, however, there might be various effects, such as data having noise and imperfect geometry which might be of importance, and which could also influence the quality of the reconstruction of real data.

Discretization of the Projection

Figure 17 shows a schematic view of the discretization of coordinates u_m und V_n in a fixed coordinate system of the body under an angle ψ_j on an element Xj of the detector 12 in the fixed coordinate system of the system.

A set of measured data will be measured for a limited number of angles ψ_j wherein j e {1 , ..., N_j}. Hereby, the set of data of each individual angle, i.e. the intensity data set generated under a specific radiation angle, comprises the grey values of a line detector or a line of a detector, having Nj elements. For simplicity and without constraint of the generality, such a vector element has a width of one unique length. The element having the index i therefore ranges from XM to Xj = -N/2+i, wherein i € (1 , ..., NJ.

Accordingly, the area having a size of N_m x N_n area units in the u-v-plane is subdivided around the origin in as much unit areas, indicated with m and n. The corresponding unit area ranges from u_m.i to u_m = -Nm/2+m and from v_n.i to V_n = - N_n/2+n.

The fundamental relation x = u cos(φ)+v sin(φ) between the coordinate system fixed in space and the coordinate system fixed relative to the body, transfers directly to the discretisized data: the simplest correspondence includes assigning to each unit area the detector element, which center is closest to the center of the area. Said correspondence is called "next neighbour correspondence", wherein:

x = (- ^y∞ + m - ύ • ∞s(φ) + (- ψ + n - \ ) • sintø) i = ψ +x+ l

and wherein i is rounded down to the next natural number. In case Nj = N_n = N_m, the relation can be simplified to

i = ^i • (1 — cos(φ) — sin(φ)) + m • cosfø) + n • sin(</?) » f + (m - N/2) • cosfø) + {n - N/2) • sinfø)

Later, a better correspondence will be described (see e.g. figure 31 ) which weights the contribution of each unit area to the respective detector element according to which ratio of the area will be projected on the respective section of the x-axis.

Tomαraphic Linear Problem

The projection can be written as

and using matrix notation with a given sinogram xι and the searched unknown f_mn as

e.. — pmn f °ιj * ij Jmn The elements of the matrix P™ describe the weights with which an f_mn is projected under the angle φ_t onto the detector element Xj. It should be proportional to the area ratio of the size of the element (u_m, V_n), which falls, by projection on the x-axis, into the segment x_M to Xj of the detector.

Such a matrix PT is exemplarily shown in figure 18. In said representation, the matrix further comprises 20 x 20 sub matrices A/ with each 20 x 20 elements

⁽A^J): =^PΓ -

Certainly, also in this case, the approximation can be used, wherein each sub area is entirely corresponding to the detector element, in which area the projected center of the (sub) area lies. Using such a next-neighbour-correspondence, the entries of the matrix P™ are either 0 or 1.

The calculation complexity and the memory complexity of the linear problem Px = s shall be considered. The projection matrix P is a (Ni • N_j) x (N_m • N_n)-matrix. For simplicity: N: = Nj « N_j « N_m « N_n. However, P is only thinly populated, which means that regarding the N⁴ entries, O(N³) are different from 0, only.

Naive methods for solving a set of linear equations, however, cannot be applied, since N is large. Particularly, N typically is on the order of N « 1000. Therefore, Gauss-Seidel-elimination would use O(N⁶) operations and O(N³) of memory, for storing the matrix. Further, P very often is not exactly invertible.

According to a preferred embodiment, the above problem is solved, using an iterative method. Here, not the matrix is a data element of main importance. Rather, advantageously, the vector, which the matrix is applied to, is most relevant. For application of a matrix on a vector, said matrix does not have to be stored - it is sufficient to program the effect of the matrix on the vector.

The following exemplary pseudo code shows a rudimentary projection and backprojection applied to a vector using 0(N²) memory. Said operations can be executed within O(N³) of time.

Pseudo code:

Projection (slice[N] [N] → sinogram [N] [N])

V i, j: sinogram [i] [j] = 0 V j: co = Cos (π • j/N) si = Sin (π * j/N)

V u:

V v: i = N/2+(u-N/2) . co+(v-N/2) . si sinogram [i] [j] += slice [u] [v]/N Backproj. (sinogram [N] [N] → slice [N] [N])

V u, v: slice [u] [v] = 0 V j: co = Cos (π * j/N) si = Sin (π * j/N) V u:

V v: i = N/2+(u-N/2) • co+(v-N/2) * si slice [u] [v] += sinogram [i] [j]/N

An algebraic method according to a preferred embodiment of the invention

The direct reconstruction equation f(u, v) = P^τFP_φ(x) leads in first approximation F « 1 to unfiltered backprojection f(u, v) = P^τP_φ(x).

Moreover, in discretized form, the backprojection P^τ is the transposed of the projection: (P^r)°_mn = P™ ■ Multiplying said backprojection to S_y = P?"f_mn a conditional equation of the unfiltered backprojection is obtained:

\^r Jm'n' ^S'J ^~ V )m'n' ^rij Jmn

This is a new linear problem:

Λ mn f fO

^Λm'ri Jmn — Jm'n¹

having the matrix

Λmn I pT\ij prnn

^Αm'n' — \^r Jm'n' ^rij

with the unfiltered backprojection as the right side:

fO _ ( pT\ij Jm'n' ^* VΛ Im'n' *»J

A is a N⁴-matrix which is not thinly populated as compared to P and P^τ which can be seen from figure 19. However, A provides a solution approach and, in particular, has properties allowing the specific method therefore:

1. The matrix A projects the real space having the coordinates (u_m, V_n) onto itself.

Since thereby images are generated, the matrix A preferably represents a digital image processing filter. The effect of A is obvious, when comparing the unfiltered backprojection with the object or from figure 19. The matrix A is a diffusion filter.

In particular, figure 19 shows schematically an image processing filter A having a subdivided matrix of (exemplary) 20 x 20 sub-matrices. Each of the sub-matrices corresponds to an image element (m, n). The exemplary 20 x 20 entries of each sub-matrix determine, how much image element (m¹, n') effects image element (m, n). Obviously, each image element is primarily identically projected. Additionally, diffusion components appear.

Moreover, digital image processing filters are often translatoric invariable, i.e. typically, the entries B™_n, exclusively depend on the differences m-m¹, and n- n¹: B™_πl = b(m-m\n -n') . In the present case, as an example, this does not apply- A cannot be reduced in such a form.

2. For inversion of digital image processing filters, for instance to correct shaking or unfocussed images, iterative methods are commonly used. For example, the conjugate gradient method can be used. However, a quadratic, symmetric and positive-definite matrix is necessary. The matrix A has said properties.

Summarized, according to a preferred embodiment of the present invention, it was found out that the reconstruction problem can be written as P^TPf = f⁰ with the unfiltered backprojection f⁰ and the matrix P^TP, which is quadratic (by projecting the R^N- x R^N" onto itself), symmetric (there is provided (P^TP)^T = P^TP) and positive definite (there is provided f^TP^TPf > 0 V x ≠ 0).

Looking at the quadratic form

F{f) = ϊ f^τP^τPf - f^xrf

and at a gradient

F'(f) = ϊ (^pτp)^τf + l P^τpf - f° which vanishes at the minimum of F. Due to the symmetry of P^TP there is provided 0 = P^TPf-f°. According to a preferred embodiment of the invention, it was found out that, advantageously, minimizing F is equivalent to solving the equation P^TPf = f°. There exist a plurality of methods for minimizing F. According to the preferred embodiment of the present invention, it was found out that, advantageously, the method of conjugate gradients is most efficient. A pseudo code for such a preferred embodiment follows:

i = 0 (initialisation) t = f (or f = 0 or ...)

(the residuum d = r (initial direction)

Onew ⁼ r r i = 1 , 2, 3 ... (iteration) q = P^TPd α = δnew/(d^τq) (optimal step size) t = f¹ + αd (step to new point/iteration step) r = r - αq (the new residuum)

O₀Id ⁼ Onew

Onew — r r β = δnew/δold d = r + βd (the new direction)

As can be seen above, input of the conjugate gradient method is the unfiltered backprojection, which can be determined by backprojecting the measured sinogram directly.

During the individual iteration steps, the method applies to projection and backprojection. This can be operatively implemented. Advantageously, according to this preferred embodiment of the invention, there is no need for storing the matrices. During the iterations, f will approach the correct solution. The convergence is certain within 0(N²) iterations, in condition that projection and backprojection correspond precisely to the real experiment and in condition that the data is free of systematical errors and noise.

However, said conditions need not be exactly fulfilled. Actually, advantageously, the method converges already after few iterations providing a very good approximation. Advantageously, according to this preferred embodiment of the invention, there is provided an iterative method/algorithm having time complexity O(N³) and memory complexity 0(H²).

For testing the new algorithm/method, the Shepp-Logan-Phantom is provided, wherein Nj = N_m = N_n = 268 and Nj = 800.

The algorithm is started with an iteration using unfiltered backprojection f°. Some of the following iteration steps f\ i.e. the reconstructed tomographic image data, are exemplary shown in figure 2OA. Accordingly, similar exemplary views are shown for a turbine blade in figures 2OB and 2OC. On the upper left hand side of figure 2OA, there is shown the unfiltered backprojection, and subsequently iteration steps , f², f⁵, f²⁰ and f¹⁰⁰. Finally, on the lower right hand side, there is shown the filtered backprojection as a comparison. Accordingly, figure 2OB shows exemplary tomographic image data of an turbine blade, wherein on the upper left hand side the backprojection is shown and on the lower right hand side the tomographic image data is shown after 20 iterations.

Figure 2OC is similar to figure 2OB, wherein figure 2OC shows a detailed section of figure 2OB.

In order to provide figures 2OB and 2OC, a resolution of Nj = N_m = N_n = 640 was used and Nj = 2000.

The quality of the image data according to the preferred embodiment of the present invention becomes obvious when directly comparing it with the filtered backprojection. Figure 21 shows the differences between the tomographic image data and the real Shepp-Logan-Phantom in falsecolours, here for simplicity shown in grey scale. The right hand side shows the difference between the real Shepp-Logan- Phantom and the reconstruction using the method according to a preferred embodiment of the invention after 100 iterations. The left hand side shows the difference between the real Shepp-Logan-Phantom and the unfiltered backprojection. The differences are indicated by means of the bar in the lower half of figure 21 , wherein the left hand side corresponds to large errors, i.e. poor reproduction quality and the right hand side corresponds to small errors, i.e. good reproduction quality. Obviously, there is a large improvement from the naive backprojection to the method according to the preferred embodiment of the present invention.

Further to that, the optical impression should be visualized by an objective error. The absolute error during the iteration step i summed over all image elements is given by

m n

wherein f is the object to be determined. In the present case the exemplary test object, i.e. the Shepp-Logan-Phantom. Worthwhile, said error is normlized by the unfiltered backprojection:

ε' = E' /E°

Said error ε^ι is shown in figure 22, which shows the decreasing error on a double logarithmic scale. Incomplete information

Advantageously, the method according to a preferred embodiment of the invention is suitable for providing tomographic image data from incomplete information. Incomplete information comprises, e.g., the case that the projection is known only for a limited angular interval. This can be the case, when the object to be examined is surrounded by tissue which should not be penetrated by radiation or which cannot be penetrated by radiation. Examples can be mammography and dental tomography. Also, metal objects might prevent the penetration of radiation, in particular x-rays. Moreover, in technical radiography, which also includes neutron-tomography, typically there exist restrictions due to the geometry of the object or strongly absorbing areas.

A specific application in neutron-tomography is used for enhancing the spatial resolution, wherein, according to a preferred embodiment, a substantially/mere flat object is tilted by only a small angle. Therefore, the object can be arranged close to the detector. In case of a divergent ray, the spatial resolution is directly enhanced.

Actually, the scintillator can be attached to the object, in particular it can be glued to the object. Therefore, resolution can be enhanced.

In case of medical tomography, radiation time and radiation doses can be minimized, advantageously speeding up tomography.

In the following, the case of a limited, not necessary continuous angular interval should be described, as a preferred example of the invention. Such a tomography is influenced by artefacts, aliasing and blurring of the object, due to the missing projection. Conventionally, reconstruction was assumed to be possible only by an angular interval of not less than 120°.

In any case, filtered backprojection is not suitable for reconstructing tomographic data from incomplete starting data, such as incomplete sinograms. In all the above described differentiations, integrals over the angular range 2π occurred. Although, in case of parallel neutron-tomography, neglecting scattering, there is provided P_φ(χ) = P_(φ+π)(-χ) , such that the entire information is already included in the first 180°.

A further reduction of the angular information, however, leads to entirely wrong representations/image data.

However, advantageously, when deriving a preferred embodiment of the present invention, no assumptions regarding the distribution of the angles are necessary. Preferably, neither have said angles to be comprised in an entire circle, nor do they have to be continuous or not even equi-angular arranged. Each projection direction corresponds to a matrix-line or matrix-column, which is not directly connected to the other ones.

Actually, advantageously, using the method according to a preferred embodiment of the present invention, tomographic image data can be reconstructed from incomplete (input) data. The reconstruction of the first 90° of the Shepp-Logan-Phantom is exemplarily shown in figures 23 and 24. In figure 23 there is shown the filtered backprojection having an incomplete data set, in particular exclusively using the first 90°. In figure 24, there is shown the tomographic image data, according to a preferred embodiment of the present invention, using the first 90°. Therefore, figures 23 and 24 can be directly compared. Advantageously, the method according to the present invention is much better suited for reconstructing incomplete data.

Advantageously, in order to reconstruct an incomplete set of data, the method can be adapted in the projection- and backprojection routine not in the iteration routine. In pseudo code, e.g. for a continuous set of angles of yet 90°, the exemplarily amendment could be:

co = Cos (π * j / N) sin = Sin (π ^* j / N) I co = Cos (π / 2 * j / N) sin = Sin (π / 2 ^* j / N)

Figure 25 shows the iteration errors, as already displayed in figure 22, for further different angular intervals, in particular for angular interval of 180°, 90°, 30° and 10°.

Further advantageously, the method according to a preferred embodiment of the present invention actually works with much smaller apertures, as can be seen in figure 26A, figure 26B and figure 26C. On the left hand side of figure 26A, there is shown image data according to a preferred embodiment of the invention, on the right hand side, image data obtained by backprojection. In figures 26B and 26C, on the right hand side, there is shown image data according to a preferred embodiment of the invention, on the left hand side image data obtained by backprojection. In particular, from said figures, it can be seen, how superior the method according to a preferred embodiment of the present invention works, in particular for small angles, compared to filtered backprojection.

The only limit is due to the fact that for even smaller angles, the line and columns of the projection matrix become more and more similar. The matrix thereby degenerates. Hence, having a less ideal conditioning of the matrix, the errors decrease slower, as can also be seen from figure 25.

Fan Beam Geometry

A preferred embodiment of the projection is the fan beam geometry, which is preferably used in the medical field. Contrary to the parallel geometry (left hand side), when applying the fan beam geometry from a radiation source, a fan beam-like radiation is provided which can be detected by a linear or curved detector 12, as shown in figure 27A. Figure 27A particularly shows a schematical view of possible arrangements of the radiation source 10 with respect to the body 14 and the detector 12. In figure 27B, there is shown a further preferred embodiment of the present invention. In particular, figure 27B shows a radiation source 10 for radiation having fan beam geometry. Further, there is shown a body 14 which is transported on a conveyor belt 15. The object may e.g. be luggage, such as a suitcase, a bag, etc. The object 14 is transported on the moving conveyor belt 15. Thereby, the object is passing the fan-beam geometry of the radiation of the radiation source 10. The penetrated radiation is detected by the detector 12, such that intensity data of the penetrated radiation can be created. Over time, different intensity data sets (referred to in figure 27B as intensity data), in particular N intensity data sets (#1 ... #N/2 ... #N) can be created. Form said intensity data sets, the tomographic image data can be obtained. As an example, in figure 27B, the position of the object 14 over time is shown for three different time steps. The time increases from the top to the bottom in figure 27B. The object 14 moves from the left side in figure 27B to the right side in figure 27B with increasing time. Therefore, the intensity data sets are generated under different angles relative to the radiation source 10, particularly without the necessity of rotation of the object 14. As an example the embodiment according to figure 27B is a schematic representation of a security system e.g. used at an airport, a train station etc.

The necessary modifications according to a preferred embodiment of the present invention are provided by a corresponding amendment of the projection between the body-fixed and the system-fixed coordinate systems. For a equi-spatial geometry having the distance L₁ between radiation source 10 and detector 12 and the distance l_2 between the radiation source 10 and the object 14, as schematically shown in figure 28, the respective section in pseudo code could be similar to

x = (m-N/2) * co + (n-N/2) ^* si y = (m-N/2) ^* si - (n-N/2) ^* co i = N/2 + L₁ ^* x / (L₂ + y) Exemplary reconstructions for a fan beam geometry are shown in figure 29. On the left hand side, there are shown the reconstruction according to a preferred embodiments of the present invention. On the right hand side there is shown reconstructions according to the filtered backprojection. The angular interval is 180° in the upper half of the figure and 90° in the lower part of the figure. In both cases, 20 iterations were used. It should be noted that for fan beam geometry, already 180° represent a limited angular interval, since the projection is not symmetric

In the above description, a continuous angular interval has been considered, e.g. starting from 0° to 10°, 30°, 90° or 180°, respectively. In the medical field and the sciences, there exist cases, wherein not the entire angular region is available. In particular, it might be necessary, to spare out individual angles or it might not be possible, that individual angles can be used. This might be the case, since outside of the re-constructional area, there might be metal inlays or such alike. In other cases, there are exclusively limited angular regions available for example by openings in a otherwise non transparent covering.

Advantageously, according the preferred embodiment of the present invention, it is not necessary that the angles are provided in a continuous interval, that they are equi-angular distributed, etc. Each available projection provides a line-column vector, which information can be used. In figure 30 an exemplary corresponding reconstruction can be seen, which is based on two angular intervals of 5° width, each. The two angular intervals are spaced apart from each other by 90°. On the left hand side there is shown the reconstruction according to a preferred embodiment of the present invention. On the right hand side there is shown the reconstruction made by a filtered backprojection.

Using 100 iterations, there is already provided a much better quality as in case of a continuous interval of 10°, as shown in figure 26A. This is the case, since the information in angular intervals having a width of 5° and rotated by 90° is more complementary than consecutive intervals of 5°, each. Due to the different line- /column vectors, the matrix is conditioned in a better way and the convergence therefore is better.

Realistic Complications

The method according to a preferred embodiment of the present invention converges reliable to the searched solution, in case two conditions are fulfilled: the simulated projection corresponds to the real projection which has to be inverted and the detected projections must be free of errors, i.e. they must be free of systematical errors such as illinearities or statistical errors such as noise.

Typically, said conditions cannot be fulfilled entirely, which can reside in a problem due to instabilities, possibly residing from the inverse and badly pronounced problem.

Preferably, the next neighbour mapping can be replaced by a correct weighting of the areas. Thereby, a projection of the areas onto the x-axis could be exactly calculated. Additionally/alternatively, each of the area elements, which shall be projected, can be covered by an equi-distant point grid, as can be seen in figure 31. Each point is projected, having a weight, wherein the weight corresponds to 1 divided by the number of points. Since said weight is applied to the entire projection, on each detector element, a natural number of points is projected. The matrix PJ"¹ advantageously can be calculated very efficiently and stored in a compressed way.

Hereby, preferably, for each j, n and m, the index of the smallest i is stored, which obtains weights of the area (m, n). Afterwards, for each element, starting with i, the number of weights is stored. When the first element is reached, which lies outside the projected area, therefore not having received any weights, 0 is stored. Advantageously, this structure can be provided by a single, preparing projection and can be stored. In the following projections and backprojections, said structure can be used, wherein only actual entries of the respective vectors are multiplied. Advantageously, there is provided an optimal time and memory use.

Preferably, the divergence of neutron radiation can be taken into account. Further preferably, the section on the x-axis, in which the weights of an area are distributed, will be widened additionally. Said widening can be preferably done depending on the distance Li - L₂ - y of said area to the detector. Advantageously, said method/algorithm works for such divergent geometries as well as for the above indicated fan beam geometry (both geometries have parallelities/similarities).

Further preferably, small numerical errors are avoided. Since said small numerical errors can increase during the iterative process. According to the invention, it was found that in particular errors in series of the trigonometric functions at specific angles are of particular interest. Therefore, preferably, a wrap-around the sine and the cosine functions should be used for covering multiplicities of 30°. In pseudo code, this could be implemented as follows:

if (6*j == 0*N) co = 1.0 else if (6*j == 1*N) co = 0.866025403784438646763723170753 else if (6*j == 2*N) co = 0.5 else if (6*j == 3^*N) co = 0.0 else if (6*j == 4^*N) co = -0.5 else if (6*j == 5*N) co = -0.866025403784438646763723170753 ... else co = Cos(π^*j/N)

Further to that, a primary constraint of the possible precision is noise in the measuring data. The consequences of noise can be easily described in a simplified manner: The tomography of an empty space provides projections, which are constantly zero. Assuming that in a single projection a single pixel comprises a small value, as a consequence of noise, said value is backprojected as a line. In order to obtain zero in the other projections, preferably, small negative values are added to the entire area. Errors can successively be increased as high-frequency disturbances.

Since an image can be assumed as a linear superposition of an undisturbed image and a disturbed image, obviously, clearly during the iterations, decreasing errors of the undisturbed image can be overwhelmed by increasing errors of the disturbed image. Therefore, there exists one iteration with a minimal overall error. For any further or less iteration, the overall error is larger than the minimal overall error. Consequently, the less noisy the measuring data are, the more iterations can be undertaken, before reaching said minimal overall error. Moreover, the less noisy the measuring data are, the smaller said minimal overall error will be.

In figure 32 there are shown reconstructions of a noisy sinogram having an exaggerated signal to noise ratio of 10. Shown are the 5^th, the 10^th and the 50^th iteration step. Obviously, in the beginning, the tomographic image data is blurry due to dominant errors (left hand side). Later on, the tomographic image data sharpens. Simultaneously, high-frequency errors increase, which are affected by noise (middle image and right hand side image).

Real Measuring Data

The method according to a preferred embodiment of the present invention was tested with real measuring data. For real data having an entire angular interval of 180° there is provided a quality, comparable to the filtered backprojection. In particular, an exemplary aluminum sample was used, as shown in figure 33. In the beginning, the entire measuring data were reconstructed, as seen on the upper part of figure 33. The filtered backprojection and the method according to a preferred embodiment of the present invention are of similar quality. Later, only half of the data, i.e. an angular interval of 90°, was reconstructed, as can be seen in figure 33, lower section. Obviously, the method according to a preferred embodiment of the present invention is superior to the filtered backprojection, as can be seen from the image of the sample having a substantially rectangular form. Simultaneously, the contrast is higher as compared to filtered backprojection.

Summarizing, according to a preferred embodiment of the present invention, advantageously, conventionally inaccessible tissue, in particular dental areas, the chin or the chest etc. can be displayed. In natural science, insights in merely inaccessible objects could be provided. In the field of neutron-tomography the spatial resolution could be increased by an order of magnitude, since it might be possible to rotate the sample only slightly and, thereby, decrease the distance to the detector.

Moreover, further preferably, the method according to a preferred embodiment of the invention could be adapted to the precise circumstances of ray and detector thereby increasing the quality of the tomographic image data.

Summarizing, one preferred embodiment of the present invention should be described as an example, by means of a flow-diagram, as shown in figures 34A,

34B and 34B. After starting the preferred method, in a first step S1 , the sinogram s is used as an input data. Therefore, even prior to step S1 , the sinogram s can be generated using radiation, such as gamma radiation and a detector. Thereby, radiation can be shined on an object or an object can be brought into the optical path of the radiation. In particular, the sinogram can be generated by rotating the object within the radiation and detecting the intensity of the attenuated radiation after penetrating the object and/or by rotation of the radiation source and the detector.

Alternatively/additionally the sinogram s can at least partly be obtained from a database, etc. In other words, the sinogram s can at least partly be generated prior to step S1. The sinograms can be generated, using a conventional device, which does not necessarily need to be connected to an apparatus suitable for carrying out a method according to a preferred variation of the present invention. In step S2, the matrix for the projection P and the matrix for the backprojection B are calculated and stored in advance.

In a following step S3, the backprojection B is obtained, in particular calculated from the sinogram s using the corresponding matrix for the backprojection B of the sinogram s.

In step S4, the solution to be obtained, i.e. vector f (e.g. relating to the above description the attenuation distribution map), is set identical to the vector b. In the consecutive steps, the solution f shall be calculated/obtained. The solution f can be the attenuation distribution map of the object regarding one slice of the object. In other words, the solution vector f can be the attenuation distribution map of the object in one plane corresponding to one radiation direction of the penetrating radiation.

In step S5, the projection of the solution vector f will be calculated/obtained and set to a temporary variable temp2.

In step S6, a temporary variable will be calculated from the backprojection on the temporary variable temp2.

In step S7, the residuum r will be calculated from the backprojection b and a temporary variable tempi .

In step S8, the direction of the iterative step d is set to the residuum and in step S9, the error is calculated by the scalar product of the vector r, i.e. the scalar product of the residuum.

In the following step S10, the iterative loop is started with setting the iteration counter i = 1. In step S11 , the temporary variable temp2 is set equal to the value of the projection of the direction of the iterative step d.

In steps S12, S13 and S14, a conventional propagation according to the conjugate gradient method is applied.

In step S15, it is determined, whether the iteration number/counter is dividable by 10. In other words, for every 10^th iteration cycle, the steps S16, S17 and S18 are carried out. For any iteration, which is not a 10^th iteration cycle, step S19 is carried out, instead. In particular, step S19 is a fast approximation of steps S16, S17 and S18. However, in order to minimize or avoid systematical errors, steps S 16, S17 and S18 have to be applied for every 10^th iteration cycle. It is also possible that steps S16, S17 or S18 are applied for every other iteration cycle, for every 5^th iteration cycle, for every 15^th iteration cycle, etc. In other words the detailed calculation according to steps S16, S17 and S18 can be carried out according to a repetition routine different to i = 10 cycles, but i can be equal to 1 , 2, 3, 4, 5, ... etc.

Following that, steps S20 to S23 are carried out, which are core steps of the conjugate gradient method, wherein the new searching direction d is determined.

Further on, in step S24, the error value of the conjugate gradient method is provided. In case said error is not smaller than a predefined value/threshold, a further iteration will be carried out. In other words, the iteration counter i will be increased by 1 in step S25 and steps S11 to S23 are carried out successively, thereby providing the error value in step 24.

In case the error value is smaller than a predefined value, in step 26, the solution f is provided.

Instead of providing a specific threshold, which the error value has to be smaller for every iteration cycle, as described above, the solution f could be displayed and the quality of the solution could be judged manually, in every iteration cycle or every other iteration cycle or every 5^th iteration cycle or every 10^th iteration cycel, etc. Thus the iteration could be aborted in case a user considers the solution f as suitable.

The above described variables are used exemplary, only. In particular, they are used for different preferred embodiments of the invention. However, the terminology of the respective variables can be substituted within each embodiment or different terminology of the variables can be used for different embodiments. For example, instead of f(u, v) the term f or μ could be used. Instead of H^τ the term B or the term P^τ or (P^τ f_m,_n, could be used, instead of H^TH the term A could be used etc. Further to that, the specific typeface applied for the variables can be substituted without specific preferences. For example f can be used instead of / , P can be used instead of P or instead of P , etc.

Figure 35 shows an exemplary system which can be comprised by an apparatus according to a preferred embodiment of the invention. In particular, it is not necessary, that the apparatus according to a preferred embodiment of the invention comprises all elements, as shown in figure 35. Rather, the apparatus can comprise only some of the described elements. Further to that, also additional/other elements can be comprised, in particular elements conventionally known in the field of image processing, computer systems, tomography etc.

With reference to figure 35, an exemplary system for implementing the invention includes a general purpose computing device in the form of a conventional computing environment 20 (e.g. personal computer), including a processing unit 22, a system memory 24, and a system bus 26, that couples various system components including the system memory 24 to the processing unit 22. The processing unit 22 may perform arithmetic, logic and/or control operations by accessing system memory 24. The system memory 24 may store information and/or instructions for use in combination with processing unit 22. The system memory 24 may include volatile and non-volatile memory, such as random access memory (RAM) 28 and read only memory (ROM) 30. A basic input/output system (BIOS) containing the basic routines that helps to transfer information between elements within the personal computer 20, such as during start-up, may be stored in ROM 30. The system bus 26 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures.

The personal computer 20 may further include a hard disk drive 32 for reading from and writing to a hard disk (not shown), and an external disk drive 34 for reading from or writing to a removable disk 36. The removable disk may be a magnetic disk for a magnetic disk driver or an optical disk such as a CD ROM for an optical disk drive. The hard disk drive 34 and external disk drive 34 are connected to the system bus 26 by a hard disk drive interface 38 and an external disk drive interface 40, respectively. The drives and their associated computer-readable media provide nonvolatile storage of computer readable instructions, data structures, program modules and other data for the personal computer 20. The data structures may include relevant data of the implementation of the method for generating at least two- dimensional image data. The relevant data may be organized in a database, for example a relational or object database.

Although the exemplary environment described herein employs a hard disk (not shown) and an external disk 42, it should be appreciated by those skilled in the art that other types of computer readable media which can store data that is accessible by a computer, such as magnetic cassettes, flash memory cards, digital video disks, random access memories, read only memories, and the like, may also be used in the exemplary operating environment.

A number of program modules may be stored on the hard disk, external disk 42, ROM 30 or RAM 28, including an operating system (not shown), one or more application programs 44, other program modules (not shown), and program data 46. The application programs may include at least a part of the functionality as detailed in the above figures.

A user may enter commands and information, as discussed below, into the personal computer 20 through input devices such as keyboard 48 and mouse 50. Other input devices (not shown) may include a microphone (or other sensors), joystick, game pad, scanner, or the like. These and other input devices may be connected to the processing unit 22 through a serial port interface 52 that is coupled to the system bus 26, or may be collected by other interfaces, such as a parallel port interface 54, game port or a universal serial bus (USB). Further, information may be printed using printer 56. The printer 56, and other parallel input/output devices may be connected to the processing unit 22 through parallel port interface 54. A monitor 58 or other type of display device is also connected to the system bus 26 via an interface, such as a video input/output 60. In addition to the monitor, computing environment 20 may include other peripheral output devices (not shown), such as speakers or other audible output.

The computing environment 20 may communicate with other electronic devices such as a computer, telephone (wired or wireless), personal digital assistant, television, or the like. To communicate, the computer environment 20 may operate in a networked environment using connections to one or more electronic devices. Figure 35 depicts the computer environment networked with remote computer 62. The remote computer 62 may be another computing environment such as a server, a router, a network PC, a peer device or other common network node, and may include many or all of the elements described above relative to the computing environment 20. The logical connections depicted in figure 35 include a local area network (LAN) 64 and a wide area network (WAN) 66. Such networking environments are commonplace in offices, enterprise-wide computer networks, intranets and the Internet.

When used in a LAN networking environment, the computing environment 20 may be connected to the LAN 64 through a network I/O 68. When used in a WAN networking environment, the computing environment 20 may include a modem 70 or other means for establishing communications over the WAN 66. The modem 70, which may be internal or external to computing environment 20, is connected to the system bus 26 via the serial port interface 52. In a networked environment, program modules depicted relative to the computing environment 20, or portions thereof, may be stored in a remote memory storage device resident on or accessible to remote computer 62. Furthermore other data relevant to the application of the insurance claim management evaluation method (described in more detail further below) may be resident on or accessible via the remote computer 62. The data may be stored for example in an object or a relation database. It will be appreciated that the network connections shown are exemplary and other means of establishing a communications link between the electronic devices may be used.

Accordingly, the invention can be implemented in digital electronic circuitry, or in computer hardware, firmware, software, or in combinations of them. The invention can be implemented as a computer program product, i.e., a computer program tangibly embodied in an information carrier, e.g., in a machine-readable storage device or in a propagated signal, for execution by, or to control the operation of, data processing apparatus, e.g., a programmable processor, a computer, or multiple computers. A computer program can be written in any form of programming language, including compiled or interpreted languages, and it can be deployed in any form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program can be deployed to be executed on one computer or on multiple computers at one site or distributed across multiple sites and interconnected by a communication network.

Method steps of the invention can be performed by one or more programmable processors executing a computer program to perform functions of the invention by operating on input data and generating output. Method steps can also be performed by, and apparatus of the invention can be implemented as, special purpose logic circuitry, e.g., an FPGA (field programmable gate array) or an ASIC (application- specific integrated circuit).

Processors suitable for the execution of a computer program include, by way of example, both general and special purpose microprocessors, and any one or more processors of any kind of digital computer. Generally, a processor will receive instructions and data from a read-only memory or a random access memory or both. The essential elements of a computer are a processor for executing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to receive data from or transfer data to, or both, one or more mass storage devices for storing data, e.g., magnetic, magneto-optical disks, or optical disks. Information carriers suitable for embodying computer program instructions and data include all forms of non-volatile memory, including by way of example semiconductor memory devices, e.g., EPROM, EEPROM, and flash memory devices; magnetic disks such as internal hard disks and removable disks; magneto-optical disks; and CD-ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in special purpose logic circuitry.

To provide for interaction with a user, the invention can be implemented on a computer having a display device such as a CRT (cathode ray tube) or LCD (liquid crystal display) monitor for displaying information to the user and a keyboard and a pointing device such as a mouse or a trackball by which the user can provide input to the computer. Other kinds of devices can be used to provide for interaction with a user as well; for example, feedback provided to the user can be any form of sensory feedback, such as visual feedback, auditory feedback, or tactile feedback; and input from the user can be received in any form, including acoustic, speech, or tactile input.

The invention can be implemented in a computing system that includes a back-end component, e.g., as a data server, or that includes a middleware component, e.g., an application server, or that includes a front-end component, e.g., a client computer having a graphical user interface or an Web browser through which a user can interact with an implementation of the invention, or any combination of such back- end, middleware, or front-end components. The components of the system can be interconnected by any form or medium of digital data communication, e.g., a communication network. Examples of communication networks include a local area network ("LAN"), a wide area network ("WAN"), and the Internet. The computing system can include clients and servers. A client and server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other.

A number of embodiments have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the invention. Accordingly, other embodiments are within the scope of the following claims.

List of Reference Numerals

10 radiation/x-ray source

12 line detector/ detector

14 object/patient

16 grid

18 arrow

20 conventional computing environment

22 processing unit

24 system memory

26 system bus

28 random access memory (RAM)

30 read only memory (ROM)

32 hard disk drive

34 external disk drive

36 removable disk

38 hard disk drive interface

40 external disk drive interface

42 external disk

44 one or more application programs

46 program data

48 keyboard

50 mouse

52 serial port interface

54 parallel port interface

56 printer

58 monitor

60 video input/output

62 remote computer

64 local area network (LAN) wide area network (WAN) network I/O a modem

Claims

Applicant: Universitat Heidelberg"Method for generating at least two-dimensional tomographic image data, apparatus, computer program product"Our Ref.: H 3102WO - ds / alClaims

1. Method for generating at least two-dimensional tomographic image data of an object from a sinogram of the object, the method having the steps:

providing a sinogram s of an object, the object having an attenuation distribution map f(u, v);

generating the attenuation distribution map f(u, v) of the object by minimizing the quadratic equation:

F(f) = Vz f^TAf - b^τf

using an iterative method, wherein

A:= H^TH

and

b:= H¹S

and wherein

H is the projection matrix of the attenuation distribution map f(u, v) onto the sinogram s and H^τ is the backprojection matrix of the sinogram s onto the attenuation distribution map f(u, v) of the object;

deducing the at least two-dimensional tomographic image data from the attenuation distribution map.

2. Method according to claim 1 , wherein the iterative method uses conjugate gradients and/or a conjugate directions and/or steepest descent.

3. Method according to one of the preceding claims, wherein the sinogram s comprises at least two intensity data sets, wherein each intensity data set is generated from radiation penetrating the object under a specific radiation direction with respect to the object and wherein for all intensity data sets generated from radiation having respective radiation directions φ_k> φι with k ≠ I, there is provided:

ψk- φi < 180°, preferably

φ_k- cpi < 120°, further preferably

φ_k- φι < 90°, most preferably

(pk - φι < 10°.

4. Method according to any one of the preceding claims, wherein for calculating the iterative method, specific trigonometric calculation operations are used, which cancel the typical rounding errors of operating microprocessors.

5. Method according to claim 4, wherein a wrap around operation is carried out on the cos(x) and the sin(x) functions of the standard libraries.

6. Method according to any one of the preceding claims, wherein the resolution of the attenuation distribution map f(u, v) is virtually enhanced when projecting the attenuation distribution map f(u, v) onto the sinogram s.

7. Method according to claim 6, wherein each point of the attenuation distribution map f(u, v) is equally divided into a subset of virtual points and each virtual point is projected onto the sinogram s, wherein the weights for all the virtual points are equal.

8. Method according to any one of the preceding claims, wherein the projection matrix H and/or the backprojection matrix H^τ are compressed projection matrices.

9. Method according to claim 8, wherein the projection matrix H comprises entries solely for elements of the sinogram s, onto which points of the density distribution map are projected.

10. Method according to any one of the preceding claims, wherein a point radiation source or a radiation source radiating parallel radiation is used.

11. Method according to one of the preceding claims, wherein the radiation is neutron radiation or positron radiation or electron radiation or proton radiation or ultrasound or a combination thereof.

12. Method according to one of the preceding claims, wherein the radiation is electromagnetic radiation having a wavelength from approximately 5 femtometres to approximately 10 nanometres.

13. Method according to one of the preceding claims, wherein the direction of the radiation relative to the object is changed by rotating the object and/or the radiation source.

14. Method according to one of the preceding claims, wherein the tomographic image date are three-dimensional data.

15. Apparatus for generating at least two-dimensional tomographic image data of an object from a sinogram s of the object, comprising:

a sinogram provision means adapted for providing a sinogram s of an object, the object having an attenuation distribution map f(u, v);

an attenuation distribution map generating means adapted for generating the attenuation distribution map f(u, v) of the object by minimizing the quadratic equation:

F(f) = ¹/₂ f^TAf - b^Tf

using an iterative method, wherein

A:= H^TH

and

b:= H^τs

and wherein

H is the projection matrix of the attenuation distribution map f(u, v) onto the sinogram s and

H^τ is the backprojection matrix of the sinogram s onto the attenuation distribution map f(u, v) of the object and

a tomographic image data deducing means for deducing the tomographic image data from the attenuation distribution map.

16. Apparatus according to claim 15, wherein the sinogram provision means comprises a radiation source (10) and a radiation detector (12).

17. Use of an apparatus according to one of the claims 15 or 16 in a security system, in particular in an airport security system.

18. Computer program product, which, when loaded in the memory of a computer and executed by the computer performs the method according to anyone of claims 1 to 14.