WO2005076038A1

WO2005076038A1 - Method and apparatus for providing high-resolution reconstruction of an observed object

Info

Publication number: WO2005076038A1
Application number: PCT/SE2005/000118
Authority: WO
Inventors: Johan Lagerros; Ozan ÖKTEM; Ulf Skoglund; Lars-Göran ÖFVERSTEDT
Original assignee: SIDEC TECHNOLOGIES AB
Current assignee: SIDEC TECHNOLOGIES AB
Priority date: 2004-02-03
Filing date: 2005-02-02
Publication date: 2005-08-18
Anticipated expiration: 2006-08-03
Also published as: SE0400212L; SE0400212D0; SE526940C2

Abstract

Fast and accurate simulated emission tomographic data, e.g. for 3D reconstruc-tion methods, can be provided by registering data related to the object using a 2D detector with known detector response, and generating simulated data using a problem defined by means of the radiative transport equation and taking into account attenuation and scattering of emitted particles from the object, and the detector response, building a model operator, which is a sequence of composi-tions of operators including: an attenuation operator modelling the attenuation, a scattering operator for modelling the scattering and a detector response operator independent of the object, for modelling the detector response wherein each operator in the composition represents only one of the three as-pects, attenuation and scattering in the radiative transport equation and detection response. - a scattering operator for modelling the scattering - a detector response operator independent of the object, for modelling the de-tector response wherein each operator in the composition represents only one of the three aspects, attenua-tion and scattering in the radiative transport equation and detection response.

Description

Method and apparatus for providing high-resolution reconstruction of an observed object

Field of the Invention.

The present invention relates to a method and an apparatus for simulating

SPECT imaging. The invention also relates to a method for high-resolution reconstruction of an observed object of the kind disclosed in the preamble of claim 11 , and an apparatus to perform the method.

Background of the invention.

In the following, the concepts of forward and inverse problem will be defined. Even though there is no formal mathematical definition of these concepts, from an applied point of view there is a natural distinction between the forward and inverse problem. A forward problem can be characterized as a problem where the goal is to determine (or predict) the outcome when the models are applied to the measured data, whereas in an inverse problem the goal is to determine (reconstruct) the causes that, using the models, yield the measured data. Often one also makes use of the term reconstruction, which in this setting refers to the process of solving an inverse problem.

Several important applications require efficient methods that can accurately and efficiently solve the inverse problem that arises. Prominent examples of such applications are within medicine. In particular it is worth mentioning radiation treatment of humans and various medical imaging techniques such as MRI (Magnetic Resonance Imaging), Computerized Tomography (CT), Positron Emission Tomography (PET), and Single Photon Emission Computed Tomography (SPECT). Various different algorithms for solving the inverse problem and the formulation of the forward problem in each of the applications listed above can be found in subject matter related textbooks. The last three medical imaging techniques mentioned above, namely CT, PET and SPECT, are all based on the same physics. The interaction between the photons and the object is in these cases described by the radiative transport equation. If we also include a model, usually called the detector response, which ac- counts for the imperfections in the detector, then we can state the relation between the measured data and the object of interest. More precisely, in the case of SPECT the forward problem is to determine the measurements that are generated when the emission and attenuation of the object under study are known together with knowledge of the detector response. The inverse problem is to de- termine the emission, (and sometimes also the attenuation) from the measured data.

Many methods for solving inverse problems are based on solving the forward problem within an iterative scheme. An example of such an iterative method is the Comet technology, as described in the International Patent Application

W097/33255, hereby incorporated by reference, (corresponding European Patent Application EP 885 430 and Swedish Patent Application 9601229-9). Since the forward problem is solved in each iterative step, we see that such iterative methods for solving the inverse problem must make use of algorithms that can efficiently and accurately deal with the corresponding forward problem. Moreover, less iteration is needed if the modelling of the forward problem is more accurate. On the other hand, solving the forward problem with a more accurate model is more time and resource consuming. Hence, iterations become costly in terms of time and memory. Thus, the desire for increased accuracy also yields more complex and time-consuming algorithms. To a certain extent, using better hardware, such as faster processors, and larger and more efficient memory, can compensate for this additional complexity.

The method and apparatus according to the invention are designed to achieve an accurate modelling of the radiative transport equation and the detector response, which yields an accurate and efficient method for solving the forward problem in SPECT. This can then be used within an iterative scheme, such as Comet, in order to solve the inverse problem in SPECT.

Iterative methods for solving the inverse problem in SPECT currently used in clinical applications utilize simplified models of the corresponding forward problem. These methods ignore the effects of one or more of the following physical phenomena: scattering, attenuation, and detector response, scattering being the phenomena that is most often ignored. In many cases, such as in low- dose SPECT, such models of the forward problems are not accurate enough..

The main reason for these simplifications in the modelling of the forward problem is efficiency, and most attempts for more accurate models yields iterative method for the inverse problem that are too slow.

An example of prior art modelling of the forward problem can be found in [1]. The object of this document is to achieve a method that is fast enough for clinical use. The forward problem, which is the problem of solving to the radiative transport equation combined with the effects of the detector response, is described by a sequence of convolution operators. However, since detector re- sponse and attenuation are considered together, the convolution kernel that appears in the forward problem becomes object dependent and difficult to estimate.

Another technical problem that arises when one deals with inverse problems in- volving the radiative transport equation for reconstruction is the amount of data that is needed. If simplified models of the forward problem are to be acceptable when solving the inverse problem, the amount and quality of the measured data needs to be rather high. However, it is often not possible or advisable to collect the necessary amount of good quality data. Typically good quality data trans- lates to a high dose of nuclear radiation matter in the object to be reconstructed, which can have a harmful effect on the object. An extreme example of this is found in electron tomography where a single high dose measurement destroys the object, which are biological macromolecules, under study. Within nuclear medicine there is also a natural desire to minimize the total exposure of the ob- ject, which now is an organ in the human body.

Other examples of restrictions in the data collection that leads to less data are time constraints and motion compensation in the data collection. The first issue is simply related to the fact that a patient cannot be fixed for more than 20 min, so a SPECT data collection should not take more than 20 min. Moreover, if the patient's heart is to be imaged in SPECT, then during the data collection one needs to compensate for the heart movements. Gating, i.e. by synchronizing the data collection with the heartbeat does this.

To summarize, in general an accurate modelling of the forward problem in emission tomography (such as SPECT and PET) enables one to use less data with lower quality while maintaining the quality of the reconstruction.

Object of the Invention

It is an object of the invention to provide fast and accurate simulated emission tomographic data having application in 3D reconstruction methods.

Summary of the Invention

This object is achieved according to the invention by a method for generating simu- lated emission tomographic data of an object, said method comprising the steps

- knowledge of the attenuation and scattering properties of the object,

- using the radiative transport equation to model the interaction between the photons and the object, said radiative transport equation taking into account the attenuation and scattering properties of the object, - using knowledge of the detector in order to define the detector response, said method comprising step of building a model operator based on the radiative transport equation, the detector response, and the experimental setting, wherein

- said model operator is a sequence of compositions of operators including: - an attenuation operator modelling the attenuation - a scattering operator for modelling the scattering - a detector response operator independent of the object, for modelling the detector response wherein - each operator in the composition represents only one of the three aspects, attenuation and scattering in the radiative transport equation, and detection response. The object is also achieved by an apparatus for generating simulated emission tomographic data of an object, said apparatus being arranged to use the radiative transport equation together with the detector response to formulate the problem to be solved in order to generate simulated data and characterized in that it comprises

- model operator means for building a model operator based on a transport equation, the detector response, and the experimental setting, said model operator means comprising - Attenuation operator means for generating an attenuation operator

- Scattering operator means for generating a scattering operator,

- Detector response operator means for modelling a detector response operator wherein each operator in the composition represents only one of the three aspects, attenuation, scattering in the radiative transport equation and detector response.

According to the invention, both scattering, attenuation and detector response are taken into account when generating the simulated data, but in a way that enables efficient handling of each of the three aspects. The method and apparatus of the invention thus provide a better model of the signals given the constraints usually imposed in emission tomographic imaging. This can be used in order to reduce noise, by improving the accuracy of the solution of the forward problem in reconstruction methods such as COMET. Thus the method according to the invention provides a method for fast solution of the forward problem when the model is given by the radiative transport equation together with the detector re- sponse.

As already mentioned one of the main applications for such a method and apparatus is within iterative methods where the goal is to solve the corresponding inverse problem. Such a method for a fast and accurate solution to the forward problem can be used in iterative reconstruction methods, such as the Comet technology, in order to create an imaging system. With the method of the invention, the effects of scattering, attenuation, and detector response are considered when solving the forward problem. The invention therefore enables iterative methods, such as Comet, to create an accurate final reconstruction of a moni- tored object within a reasonable timeframe.

According to a preferred embodiment the apparatus is arranged to work with the radiative transport equation that has been simplified by the following assumptions:

- that each particle is scattered a limited number of times, for example, at most once (single scattering),

- that the scattering angle is small,

- that only particles that reach the detector in a direction substantially parallel to the detector normal will be detected (narrow beam).

These simplifications reduce the computational load without affecting the quality of the result significantly for emission tomographic imaging.

Preferably, the attenuation operator is generated based on an attenuation map, the scattering operator is generated based on the attenuation map and the detector re- sponse operator is generated based on the properties of the detector. The apparatus can further comprise attenuation map generating means arranged to receive transmission data from the camera system and to use these transmission data to generate the attenuation map.

The attenuation operator means may be arranged to receive said attenuation map from the attenuation map generating means and to generate said attenuation operator from said attenuation map.

In one embodiment, the attenuation operator is generated by performing a cone beam transform on the attenuation map.

The attenuation map may be based on a reconstruction obtained from a transmission experiment, or on experience and/or prior knowledge.

The scattering operator may be defined as a 3D reconstruction with a kernel derived from the attenuation map.

The detector response may be defined as a 3D convolution evaluated in a 2D plane with a distance independent 3D-kernel.

Alternatively, the detector response may be defined as a 2D convolution in a 2D plane, with a distance dependent 2D kernel, further comprising the step of - making a virtual slicing of the object parallel to the detector plane - assigning a distance to the detector for the slice - treating the detector response for the virtual slice as a 2D convolution with a convolution kernel dependent on the distance to and direction of the 2D de- tector - combining the treated detector responses into a 3D convolution evaluated on the detector plane.

In this case, the apparatus comprises - means for making a virtual slicing of the object parallel to the detector plane - means for assigning a distance to the detector for the slice - means for treating the detector response for the virtual slice as a 2D convolution with a convolution kernel dependent on the distance and direction to the 2D detector - means for combining the treated detector responses a 3D convolution evaluated on the detector plane. The effects of scattering, attenuation, and detector response in the solution of the radiative transport equation in 3D are described by direction dependent sequence of 3D-convolution operators. More precisely, for a fixed direction (which is the direction of the detector), the function representing the 3D object is "weighted" with the attenuation map. This requires an efficient method, e.g. ray tracing or convolution techniques, for calculating the cone-beam transform in 3D of the attenuation map.

Once the attenuation map is included for a fixed direction into the function representing the 3D object, one must also compensate for the effects of scattering and detector response. These can be modelled as 3D-convolutions with suitable kernels. First one applies the 3D convolution representing the scattering where the kernel is the scattering kernel. Then, the resulting function (which is still a function in 3D) is convolved in 3D with kernel representing the detector response. This latter 3D convolution is only evaluated for points in the detector plane. Moreover, in most models for the detector response the convolution with the detector response kernel splits into two parts, the first, which is a true 3D- convolution, the second, which is a pure line integral calculation, i.e. the value of the X-ray transform in 3D. The X-ray transform can be evaluated by well- known methods. Many similar methods, all based on some form of ray tracing and/or volume splatting, are reviewed, for example, in [2] for n = 2 and [3] for n = 3 . A different approach based on decomposing an arbitrary line into a linear combination of short line segments is given in [4] for n = 2 and partly extended to n = 3 in [5]. Another popular method, see e.g. [6] and [7], is based on the Fourier slice theorem for the X-ray transform. This requires an interpolation in Fourier space, which is known to be problematic. For some recent developments reference is made to [8].

The above method can be used according to the invention as part of a method for high-resolution reconstruction of an object in which an iterative reconstruction method is used for the reconstruction, comprising the steps of

- detecting in at least one 2D detector measured data about particles emitted from the object - creating a model operator by performing the method as defined above,

- supplying said model operator to a device arranged to calculate a goodness of fit measure between the measured data and synthetic data generated by the model operator.

This method for high-resolution reconstruction can further comprise the step of generating a first image on the basis of said information about particles.

Similarly, the above apparatus can be used according to the invention in or together with an apparatus for high-resolution reconstruction of an object, arranged to use an iterative reconstruction method for the reconstruction, said apparatus comprising: - camera means for receiving from at least one 2D detector information about particles emitted from the object, - an apparatus according to the above, arranged to receive camera parameters and transmission data from said camera means, and - model operator means for supplying said model operator to a device arranged to calculate a goodness of fit measure between the measured data and synthetic data generated by the model operator.

This apparatus can further comprise filtered backprojection means for generating a first image on the basis of said transmission data received from the camera means. Digital image processing and image reconstruction, using the inventive methodology enables one to either obtain more accurate final reconstruction of the object of interest, or receive substantially the same quality in the reconstruction using a smaller amount of measured data.

The inventive method and system can be used for recording images of living tissue, so called gated SPECT, providing a higher image quality than earlier using the same amount of radioactive nuclides as necessary for prior art 3D-imaging. Correspond- ingly, the inventive method and system can be used for recording images of living tissue for receiving substantially the same image quality, however with less noise, using a smaller amount of radioactive nuclides than is necessary for prior art 3D- imaging.

References

[1] Beekman et al., Object shape dependent PSF model for SPECT imaging, IEEE Transactions in Nuclear Science, 401, pages 31 — 39, 1993. [2] Zhuang et al., Numerical evaluation of methods for computing tomographic projections, IEEE Transactions on Nuclear Science, 41 :4, pages 1660 — 1665, 1994.

[3] Guerrero et al., Evaluation and application ofreprojection methods for 3D PET, IEEE Conference Record Nuclear Science Symposium and Medical Imaging Conference, 204, pages 1101—1105, 1993. [4] Brandt et al., Fast calculation of multiple line integrals, SIAM Journal of Scientific Computing, 204, pages 1317—1429, 1999.

[5] Boag et al., A multilevel domain decomposition algorithm for fast 0(N² logN) reprojection of tomographic images, IEEE Transactions of Image processing, 99, pages 1573—1582, 2000. [6] Westenberg et al., Frequency domain volume rendering by the wavelet X-ray transform, IEEE Transactions on Image Processing, 97, pages 1249 — 1261, 2000.

[7] Westenberg et al., An extension of Fourier wavelet volume rendering by view interpolation, Journal of Mathematical Imaging and Vision, 14, pages 103—1 15, 2001.

[8] Fourmont, Non-equispaced fast Fourier transform with applications to tomography, The Journal of Fourier Analysis and Applications, 9:5, pages 431 — 450, 2003.

Brief Description of the Drawings

For a more complete understanding of the present invention reference is now made to the following description of examples of embodiments thereof- as shown in the accompanying drawings, in which: Fig. 1 shows briefly a first embodiment of a system, where the invention can be used;

Fig. 2 shows briefly a block schedule of a preferred embodiment of an apparatus for data acquisition and processing according to the present invention, used in order to achieve a 3D-image reconstruction of an object;

Fig. 3 shows a part of the block schedule in Fig. 2, in greater detail, related to the invention as applied to SPECT;

Detailed description of preferred embodiments

Fig 1 illustrates an apparatus for providing SPECT data, known per se. In most cases one is interested in a subset H, also called the region of interest, of the object P, here shown as an elliptic shape. In medical SPECT imaging the object P is a human and the region of interest H represents some organ, e.g. the heart. In SPECT, the object P gets an injection by some gamma-radiant drug, which es- pecially is absorbed by the region of interest H. The object P is then kept fixed in a position such that the gamma camera can record data emitted from it, and especially from the region of interest H.

Typically, detectors 1 and 2 are fixedly mounted on a moving setup 9. The setup 9 can be rotated to wanted positions. It can also be moved in directions normal to the plane of rotation. The size of the detectors must be such that they are able to measure all emission data from the region of interest H at a fixed position. Preferably it must be able to measure all emission data from the object P.

Detectors 1 and 2 are each made up of a collimator 4 and a gamma camera 3.

The collimator 4 comprises cylindrical pipes arranged as a square or rectangular matrix in a substantially perpendicular direction with respect to the gamma camera 3. They are preferably manufactured in lead, Pb, in order to be used as "lenses" for gamma-photons 5 incoming in a direction normal to the gamma camera 3. A gamma photon that reaches the gamma camera is absorbed by the crystal, producing a light flash, which is detected by photo-multiplier tubes 6. When these light photons are detected, an electrical signal is generated, representing the 2D-position in the crystal of the incoming gamma photon derived as a mean of the detected light photons.

This signal is further processed in a preliminary processing means 7. The processed signal from the means 7 is then fed to and further processed in a system processing means 8, such as a computer. The computed signal is registered as a 2D coordinate c _k where j is the current position of the detector 1 or 2 itself and k indicates the 2D pixel in the detector 1 or 2.

Additional instrumentation can be included in order to measure gated data, which is of essence when the region of interest H is a fast regularly moving organ such as the heart. In order to make a recording of the heart movement the heartbeat is divided into several, e.g. eight, phases. In addition, at least one sensor E is placed on the pa- tient's body in order to register the electrocardiogram of the heartbeats. The skilled person knows how to register gated data.

The camera system 10 should also be provided with some simple camera (not shown) to provide a prior prejudice distribution. This simple camera means provides a blurred image of the gamma-radiation from the object H in a conventional way, which does not give a quantitatively correct reconstruction of the measured object but gives a qualified initial guess of the distribution of the object measures. It is to be noted that for gated SPECT a different prejudice distribution is recorded for each phase of a heart beat.

The different 3D-images recorded from the moving heart are processed in the same way. Thus, the recordings for each of the 3D-images are grouped together and then each of the images grouped in this way is processed separately.

The gamma camera system 10 is adapted for communication with the control computer 8, to which sample data are sent. The sample data comprise information about the photons detected in the photo multipliers, and are processed from tabulated form to histograms for different time intervals, camera angles and focal planes.

The control computer 8 is a part of or is adapted for communication with the computer 11 for reconstruction; to which said histograms are transferred. The camera system 10 delivers as outputs emission data, camera parameters, and transmission data.

Delivered information is herein processed as to achieve a 3D reconstruction of the observed sample object, in accordance with the present invention. The 3D reconstruction may be further processed to desired format, such as so-called Bulls eye plots, and showed by presentation means 30. As mentioned above the radiative transport equation is essential to the invention. The radiative transport equation as such is well-known physics, described in advanced physical textbooks, but its application according to the invention will now be made apparent in relation to the implementation embodiment shown in Figures 2 and 3.

The radiative transport equation.

In this document, the symbols 9?, ^ςR⁺, and 9ϊ³ denote the set of real numbers, the set of positive real numbers, and the Euclidian three-dimensional real vector space over R, respectively. Also, let S denote the unit sphere in SR³.

In the following, the building blocks of the radiative transport equation will be discussed. Given a compact (^-manifold rc5R³, let E_r(t) denote the total energy flux through r at time t. Then for each fix t the intensity is a function I, : S² x 9ϊ³ — »9T implicitly defined through the relation

where, n(T,x) denotes the unit normal of T at x and dω is the measure on T. The physical interpretation of I_t( o,x) is that it measures the intensity of the energy/particles at x in the direction ω at time t. Note that intensity I_t depends on the point x and the direction ω, and the pair (ω,x) encodes a ray (half-line) in 9?³ that ends in x and has direction ω. Since we will keep the time variable t fix, we will henceforth suppress the dependency on t.

For a fix point x e $R³ and a direction ω e S², j( o,x) describes the total emission at x in the direction ω, and is given as

j(ω,x) = ;_emission (ω,x) + ;^' _scatter (ω,x) (2) where (co,x) represents the emission of particles at x in the direction ω and j_∞Ma(ω,χ) represents the scattering of particles at x into the direction ω. The latter function can be expressed as _ca„_er (*>, ) = f_π ^~ \_s, I(ω,x)Φ(ω,ω x)dω^< (3)

where σ is the scattering coefficient, and Φ is the scattering phase function that is normalized as _s2 Φ(ω,ω x)dω'= \. (A)

The scattering phase function (ω,ω x) represents the ratio of the particles that scatters at x from a beam with direction ώ into a beam with direction ω. The motivation behind the above relation is that the intensity of particles scattered at x into direction ω equals the sum/integral over all directions of e S² of intensities I(ω',x) of particles from direction d that scatters at x into direction ω. One must therefore weigh each intensity I(ω',x) in this sum/integral with the ration <£(ω,ύ ,x) of the particles that travel in direction d and scatter at x into direction ω . The material is said to be forward scattering at x, if (ca,aj,x) «0 whenever the angle between the directions ω and ώ is larger than some small fix angle. Finally, we let β(ω,x) > 0 denote the extinction coefficient at x in the direction co and it equals

β = κ+ σ (5)

where κ(ω,x) denotes the absorption coefficient at x in the direction ω.

We now set up the radiative transport equation at x in the direction ω with ini- tial intensity I₀(ω,x) :

I ω ■ V (ω,x) = -β(ω,x)I(ω,x) + j(ω,x) 1 I₀(ω,x) := lim_i→∞ I(ω,x + sω) The factor I₀(ω,x) involved in the boundary condition of the above integro- differential equation represents the intensity given at the "beginning" of the beam (ω,χ), which is the beam that ends in x and has direction ω. Also, observe that

ω ■ V (ω,x) = — (ω,x + sω) (7) ds s=0

We will assume that the emission and the extinction are isotropic, i.e. the functions j_emssι0n and β do not depend on the direction ω. Then we define

/(*) -Λπ.ss.on *) (^keeP ^{in mind that} _mιss,on (*>> ) does not depend on ω) and we also suppress the dependency of β in ω. In this setting the radiative trans- port equation (3) can be written as

A more accurate model of the measurement process, which also includes a model of the effects of the detector, leads to the definition below.

T(f,β)(n,x) := \_sl η(ω,n,x)l(ω,x)dω (9)

where η is the detector response function and η(ω,n,x) measures the probability that a particle arriving at x from direction ω is detected when the detector has direction n at x. Note that x is a point on the detector, since we are only interested of the solution to the radiative transport equation for the area of the detector.

Finally, let us discuss the physical units of the above quantities. We begin by presenting the units in question. First, we have the unit J (Joule) for energy. Observe that this unit can be translated/interpreted as number of particles if the particles all have the "same" energy and the same charge. The unit for the solid an- gle is called steradian and is denoted by s_r. The solid angle corresponding to all of three-dimensional space being subtended is Aπ steradians and the operation of integration over S² adds the corresponding unit s_r to the expression. Finally, m (meters) denotes the unit for length, and s (seconds) is the unit for time. This gives the following table of physical units for the components of the transport equation:

Symbol Unit E_r(t) J s^"1 I,(ω,x) J s m ² s_r ^_1 Table 1. Φ(ω,ω',x) -1 Table of J^'(ω,x) J s^'1 s_r ^"' m^-3 the units σ(x),β(ω,x), κ(ω,x) m for the components of the transport equation.

The inverse and forward problems.

The forward problem for equation (8) is to calculate T(f,β)(n,x) for a given point (n,x) when and β are known. The inverse problem for equation (8) is to recover either / or β, or both given information about T(f,β) . Assuming that equation (8) can be solved for each (n,x), the solution T(f,β)(n,x) will be an expression containing / and β. Thus, the inverse problem is to recover either / or β, or both in 9t³ given the measurements

C_{j k} := T(f,β)(n_j,x_k) + measurement error (10)

for a finite set (n x_k) e S χ 9r with j = l,...,N and k = l,. m.

In CT and ET there is no emission, i.e. ≡ 0 and the inverse problem of interest is to calculate β given the measurements c_{J k}. In PET and SPECT the inverse problem of interest is to calculate / given the attenuation β and the measurements C , .

Convolution form of the solution of the radiative transport equation.

A central problem of the present invention is the forward projection problem for equation (8) and how it enables solving the forward problem in an efficient manner. The underlying basic idea is, under some mild simplifying assumptions, to express the operator T inequation (9) in terms of convolution operators.

A fast solution of the forward problem can be found by simplifying the radiative transport equation (8) so that one can solve it and express the solution / by means of convolutions.

To enable the fast solution of the forward problem some minor simplifications are made: First, it is assumed that the phase function Φ can be written as

Φ(ω,ω',x) = φ(ω,ω'). (1 1)

This simply means that the ratio Φ(ω,ω',x) of particles with direction ω' that scatter at x into direction ω can be divided into a part φ(ω,ω') independent on the point x. It is also assumed that the detector response function η does not depend on x, which is a point on the detector. Therefore, the x variable can be suppressed in η. For further simplification it is assumed that I₀(co,x) = 0, that it, that there is no outside radiation.

As already mentioned, the above simplifications are of minor nature. We will make three more assumptions that are important to the invention. These assumptions are as follows:

Single scattering: This means that particles are scattered at most one time. This assumption is justified since the probability of a particle being scattered two or more times is very low.

Forward scattering:

Forward scattering implies that φ(ω,ω') » 0 whenever ω deviates from ώ¹ by more than a small angle. Equivalently, when scattering does occur from direction ω' into direction ω, then we can assume that the angle between ω and ω' is small. Any photon that is scattered at a too large angle will be filtered out be- cause its energy is too low. The size of the angles for photons that are not filtered out are therefore related to the setting of the energy filter. For most common settings, an angular range of at most 10 degrees, i.e. an angle of 5 degrees, is considered small.

Narrow beam:

The narrow beam assumption says that only particles that reach the detector in direction almost parallel to the detector normal will be detected. More precisely, if n denotes the normal to the detector plane, then only particles that reach the detector with direction ω where ω and -n have a small angle is detected.

It turns out that it is more practical to emphasize the dependence of / in ' s_sc_ι atter > which is given inequation (9), so we change notation by replacing y_scatter with H(I). If is a fixed point and particles that reach x with the direction ω' are considered, then only two things can happen. Either a particle scatters at x or it does not. If it scatters at x then, due to the scattering assumption, the particle has not scattered before reaching x and will not scatter after passing x . Hence the scattering term H(T)(ω,x) in the radiative transport equation (8) at (ω,x) can be written as j_scma(ω,x) = H(f)(ω,x) where H is given by (9) and f is the solution to the radiative transport equation without scattering. Thus, ∞. - β(x-tω)dt f(ω,x) = ) f(x + sω)e ° ds (12)

and the radiative transport equation (8) becomes

which can be solved as I(ω,x) = I^* (ω,x) + R(ω,x) where R is defined as

∞. -j β{x-ιω)dt R(ω,x) := } H(f)(ω,x - ω)e ^• ds (14)

Inserting the definitions of into H(f) and the single scattering, forward scattering, and narrow beam assumptions, justify the following approximation:

R(ω,x) w -— J σ(x + sω){f_atlemated(-n,-) * k(-n,-ftx - ω)ds (15) Aπ ₀

where

(16)

Using the above simplification for R(ω,x), the solution I(ω,x) of equation (13) can now be expressed as

I(ω,x) ∞ f(ω,x) + ~ 1 7 ) σ(x - sω)(f_Memateά(-n,-) * k(-n,-))(x -sω)ds (17) Aπ ₀

where ω is close to -n which is the normal to the detector plane.

If one has an ideal detector, then the operator T that yields the solution to the forward problem is given as T(f,β)(n,x) = I(-n,x) . But we must include the detector response, which is given as a convolution with the function η , and we as- sume that η does not depend on the point x , i.e. the detector behaves in the same way all over the detector. Hence by equation (9)

T(f,β)(n,x):= _ιη(ω,n)I(ω,x)dω. (18)

Again, inserting the definition of /' into the above expression and utilizing our simplifying assumptions and using the fact that x is outside the support of the object / for which data should be simulated and β allows us to make the following approximation: _slf(ω,x)η(ω, )dω~\_sl \f(x-sώ)e- η(ω,n)dsdω(x) ^

^: (/attenuated (~> ) → ))

where κ(ω,y)

Moreover, ₂ - sω){f_Menmted (-«,-) *

- sω)η(ω, n)dsdω =

{σ( (/a„eπated (^~> ) * *(-* ) *

Summarizing, for each point (n,x) e S² x 5H where n denotes the normal to the detector plane and x is a point on this two-dimensional detector plane,

where /attenuated fø j ^•"=

k(ω,y) -~ φ(ω,yl\y\]y\^~2 κ(ω,y) := η(y/\y\,ω]y\^'2.

This simplified radiative transport equation (21) can be considered for each point on the detector. As can be seen, the equation comprises a series of two convolutions, each denoted by *. Such convolutions can be done in a fast and efficient manner, using hardware- or software-based technology.

The model operator is found by applying T in each of the points (n,x) that are given by the experimental setting.

The detector response has not yet been expressed in a feasible way. To enable a solution this must be done in the following.

The detector response

A common approach, see e.g. [1], for dealing with the detector response is to "slice" the object in an imagined procedure, such that each thought "slice" has a different distance to the detector. In this way one can obtain a simplified model for incorporating the detector response as a 2D convolution with a kernel that depends on the detector direction and distance to the slice. The convolution for the detector response could then be treated separately for each "slice" as a 2D- convolution with this kernel, which depends on the direction and the distance between the detector and the actual "slice". However, if this method is to be accurate, the number of slices must be large which yields a method that is not efficient enough for creating 3D-images in reasonable time.

A more efficient method, and in accordance with the invention, is to model the detector response as above, i.e. as a convolution in the 3D-space, however with values evaluated at the 2D-detector-plane. This corresponds to taking an infinite number of infinitely thin "slices" in the 2D-detector response model mentioned above.

As we have seen, the convolution kernel K for the detector response can be de- composed into one part relating to the distance between the detector and the object and one part that only deals with the detector itself, η :

where ω gives the direction of the detector in 9?³ and the pair (ω,y) encodes a ray in 9l³ that ends in the point y , which is a 2D-point in the detector plane, with direction ω . The |-|^" part of the kernel K relates to a distance between the detector and the object and the η function relates to the detector itself. This function is from now on referred to as the detector response and models the detector itself, which in the PET and SPECT case is the effects of the collimator and the gamma-camera.

If we wish to utilize the detector response kernels from such two-dimensional models, we need to find the relation between our 3D-convolution kernel κ(ω,-) and the 2D-convolution kernel that we wish to use. In general, such a 2D- convolution kernel is given as ψ : 91 x 91² → 91 , where the 91 -variable represents the distance to the detector and 91² -variable a point in the detector plane. To express the relation between these two kernels, we first embed the 9?² -variable of ψ into 91³ , i.e. define ^ : 9? χ 9?³ -> 91 as

φ(s,x) = ψ(s,x') (23)

where x e 9?³ is a point in the detector plane and x'β 9?² is the corresponding two-dimensional coordinates in the detector plane of the same point. Then, k(ω, y) = φ{y - ω,y - (y - ω)ω). (24)

In most models for the detector response the convolution with the detector response kernel splits into two parts, the first which is a true 3D-convolution and the second which is a pure line integral, i.e. the value of the X-ray transform in 3D. The true 3D-convolution can be calculated by Fourier techniques. For the second part, which is the pure line integral, there are many different algorithms.

Application in an Imaging System

Figure 2 is a logical diagram of the iterative reconstruction method (inverse problem solver) Comet for solving the inverse problem in emission tomography. It also incorporates the invention which is used in order to solve the forward problem. Although the discussion, to some extent refers to Figure 1, which describes SPECT imaging in particular, the skilled person will understand how to apply the diagram of Figure 2 to other tomographic inverse problems, such as PET, CT, and ET (Electron Tomography).

Figure 2 shows the camera system 10 of Figure 1, adapted for communication with the presentation means 30 through a number of logical blocks. It is to be noted that different steps in a computer program could provide the different measures shown in Figure 2. A dot-dashed line I-I separates the part of the diagram of Figure 2 that belong to the invention from the ones that are prior art. The inventive blocks are found to the left of the line I-I, whereas the prior art imaging system, apart from the camera system is found to the right of the line I-I. The prior art system is just an example, and could be replaced by any prior art imaging system. The circuitry to the right of the line I-I in Figure 2 is mainly based on the teachings of EP-0885439, which hereby is incoφorated by reference. The computer provided with the stated program could be a part of the reconstruction computer belonging to the camera system 10. In the example discussed in connection with Figure 1, of imaging a moving object, such as the heart, the circuitry in Figure 2 could be made in the same number as the phases of the heart beat to be shown as a separate 3D-image. However, it is also possible, and preferable, to make the different 3D-images in a sandwiched process using the same circuitry but synchronized with the phases of the heartbeat.

According to the invention, two new blocks are introduced: an attenuation map block 40 and a model block 50, both of which are connected to the camera system 10. The attenuation map may be generated based on a reconstruction ob- tained from a transmission experiment and provided to the attenuation map block from the camera system 10, or on experience and/or prior knowledge of the object.

The attenuation map block 40 is arranged to receive transmission data from the camera system 10 and to use these transmission data to generate an attenuation map, as discussed above. For a fixed direction (which is the direction of the detector), the function representing the 3D object is "weighted" with the attenuation map, generated in a block 40. The model block 50 is arranged to receive camera parameters from the camera system 10 and also the attenuation map from the attenuation map block 40, and generate a model of experiment based on the parameters and the attenuation map. The effects can there be modelled as 3D-convolutions with suitable kernels.

The generation of the attenuation map requires an efficient method, e.g. ray tracing or convolution techniques, for calculating the cone-beam transform in 3D of the attenuation map. Once the attenuation map is included into the function representing the 3D object, one must also compensate for the effects of scattering and detector response. In Figure 2 the camera system 10 is connected to three blocks that are arranged to receive emission data from the camera system 10. These blocks are all known per se. An initial guess block 60 is arranged to receive emission data from the camera system 10, which can be used as the prior prejudice distribution. The initial distribu- tion output from the initial guess block 60 is fed to a distribution block 92 which is arranged to assume a distribution which, at the beginning of an iteration process is represented by the initial guess. The assumed distribution from the distribution block 92 is fed to a first input of a first computation block 91.

A prior distribution block 70, is arranged to receive emission data from the camera system 10, and generate digitized data of an estimated reconstruction, e.g. a 3D reconstruction made with conventional methods, which is used as the starting value for the iterations. The prior distribution, generated by the prior distribution block 70, is fed to an information measure block 95, which generates an information measure operator, which is fed a first input of a second computation block 96. The distribution output of the distribution block 92 is fed to a second input of the computation block 96. The second computation block 96 is arranged to compute an information measure, which is fed to a first input of the quality of distribution block 94.

Further, an estimate block 80, is arranged to receive emission data from the camera system 10 and generate a variance estimate for individual observation grid points in each recorded observed data for the object to be imagined. The variance estimate is fed to a goodness of fit block 90, arranged to receive data also from the model block 50. The goodness of fit block 90 generates a goodness of fit operator, which is fed to a second input of the first computation block 91. The first computation block 91 is arranged to compute a goodness of fit, i.e. how near to the measuring data the output from the goodness of fit block 90 is. For this, it uses the information on its two inputs, i.e. the data received from the distribution block 92 and the goodness of fit block 90. The goodness of fit computed by the second computation block 91 is fed to a second input of a quality of distribution block 94, which determines the quality of the distribution.

The quality of distribution block 94 determines the quality of the distribution. If the quality of distribution is not good enough, the quality of distribution block 94 sends that message to an improved distribution block 97, which is arranged to derive an improved new distribution using the output from the quality of distribution block 94. This new distribution is fed as a new assumed distribution to the distribution block 92. An iteration step is provided involving the first computation block 91 and the second computation block 96. If the quality of distribution block 94 finds that the quality of distribution still is not good enough a new iteration step is provided. However, if the quality of distribution block 94 determines that the quality is adequate, it is arranged to feed a final distribution to the presentation means 30. The 3D image can then be displayed at the presentation means 30.

Figure 3 shows a more detailed embodiment of the invention, adapted to SPECT. As can be seen in Figure 3, the transmission data from the camera system 10 is first fed to a filtered backprojection (FBP) block 100, which provides a first guess in accordance with known technique. The FBP block 100 gives an image of the detection provided in relation to the object. The image generated in the FBP block is forwarded to a cone beam transform block 110 and a scattering kernel block 130.

The effects of scattering, attenuation, and detector response in the solution of the radiative transport equation in 3D are described by direction dependent se- quence of 3D-convolution operators. It is to be noted that it is only the forward scattering, i.e. towards the sensor (or sensors) 1, 2 that is to be detected. Scattering in all other directions will be ignored.

The cone beam transform block 110 is arranged to generate the attenuation map and forward it to an attenuation operator block 120 arranged to generate an attenuation operator, which is fed to a first input of a first convolution operator block 140. The scattering kernel block 130 is arranged to generate a scattering kernel, which is fed to a second input of the first convolution operator block 140. The scattering kernel block (130) is arranged to define the scattering operator as a 3D reconstruction with a kernel derived from the attenuation map. The convolution operator block 140 is arranged to generate a 3D convolution, representing the scattering where the kernel is the scattering kernel generated in the scattering kernel block 130. The 3D convolution is forwarded to a combination block 150, arranged to generate a combined attenuation and scattering operator. This attenuation and scattering operator will be fed to a first input of a second convolution operator block 170.

A PSF kernel block 160 is arranged to receive camera parameters from the camera system 1 10 and to generate a PSF kernel based on these parameters. This PSF kernel is fed to a second input of the second convolution operator block 170. In other words, the PSF kernel block is arranged to make a virtual slicing of the object parallel to the detector plane, assign a distance to the detector for the slice and treat the detector response for the virtual slice as a 2D convolution with a convolution kernel dependent on the distance and direction to the 2D detector

Then, the resulting function (which is still a function in 3D) is convolved in 3D with kernel representing the detector response. This latter 3D convolution generated in the second convolution operator block 170 is only evaluated for points in the detector plane. Thus, according to the present invention, in a block 170 using the outputs from the blocks 150 and 160, the detector response is modelled as a convolution in 3D-space, but using values evaluated at the 2D- detector plane. The block 170 generates a model operator, which is fed to the block 90, which generates goodness or fit measure also using the variances provided by the block 80 in Figure 2.

Figures 2 and 3 are logical representations. Some of the functions could be implemented in hardware, for example the block 170 for generating a model operator.

In most models for the detector response the convolution with the detector re- sponse kernel splits into two parts, the first, which is a true 3D-convolution, the second, which is a pure line integral calculation, i.e. the value of the X-ray transform in 3D. The X-ray transform can be evaluated by well-known projection techniques such as ray tracing or volume splatting or by Fourier methods using the Fourier slice theorem for the X-ray transform.

The digital image processing and image reconstruction using the methodology described above enables one to either obtain a more accurate final reconstruction of the object of interest, or receive to substantially the same quality in the reconstruction using a less amount of measured data.

It is to be noted that the gamma radiation injected in the object is mostly attenuated by the tissue, if the object is a living body, however it is also scattered in some extent.

Claims

1. A method for generating simulated emission tomographic data of an object (f), said method comprising the steps of

- obtaining or estimating knowledge of the attenuation and scattering properties of the object,

- using a radiative transport equation to model the interaction between the photons and the object, said radiative transport equation taking into account the attenuation and scattering properties of the object,

- using detector knowledge of a detector in order to define a detector response, said method comprising step of building a model operator based on the radiative transport equation, the detector response, and the experimental setting, wherein

- said model operator is a sequence of compositions of operators including: - an attenuation operator modelling the attenuation - a scattering operator for modelling the scattering - a detector response operator independent of the object, for modelling the detector response wherein

- each operator in the composition represents only one of the three aspects, attenuation and scattering in the radiative transport equation, and detection re- sponse.

2. A method according to claim 1, wherein the radiative transport equation has been simplified by the following assumptions:

- that each particle is scattered a limited number of times - that the scattering angle is small,

3. A method according to claim 2, wherein it is assumed that each particle is scat- tered at most once (single scattering).

4. A method according to any one of the preceding claims, further comprising the steps of

- generating the attenuation operator based on an attenuation map - generating the scattering operator based on the attenuation map

- generating the detector response operator based on the properties of the detector.

5. A method according to claim 4, wherein the attenuation operator is generated by performing a cone beam transform on the attenuation map.

6. A method according to claim 5, wherein the attenuation map is based on a reconstruction obtained from a transmission experiment.

7. A method according to claim 5, wherein the attenuation map is based on experi- ence and/or prior knowledge.

8. A method according to any one of the claims 6-7, wherein the scattering operator is defined as a 3D reconstruction with a kernel derived from the attenuation map.

9. A method according to claim 4, wherein the detector response is defined as a 3D convolution evaluated in a 2D plane with a distance independent 3D-kernel.

10. A method according to claim 8, wherein the detector response is defined as a 2D convolution in a 2D plane, with a distance dependent 2D kernel, further compris- ing the step of - making a virtual slicing of the object parallel to the detector plane - assigning a distance to the detector for the slice - treating the detector response for the virtual slice as a 2D convolution with a convolution kernel dependent on the distance and direction to the 2D detector - combining the treated detector responses into a 3D convolution evaluated on the detector plane.

11. A method for high-resolution reconstruction of an object in which an iterative reconstruction method is used for the reconstruction, comprising the steps of

- detecting in at least one 2D detector measured data about particles emitted from the object

- creating a model operator by performing the method according to any one of the claims 1-6, - supplying said model operator to a device arranged to calculate a goodness of fit measure between the measured data and synthetic data generated by the model operator.

12. A method according to claim 11, further comprising the step of generating a first image on the basis of said information about particles.

13. An apparatus for generating simulated emission tomographic data of an object (f), said apparatus being arranged to use the radiative transport equation combined with the detector response to formulate the problem to be solved in order to generate simulated data and characterized in that it comprises

- model operator means (100, 110, 120, 130, 140, 150, 160, 170) for building a model operator based on the radiative transport equation, the detector response, and the experimental setting, said model operator means comprising

- Attenuation operator means (110, 120) for generating an attenuation operator - Scattering operator means (130) for generating a scattering operator,

- Detector response operator means (160) for modelling a detector response operator wherein each operator in the composition represents only one of the three aspects, attenuation and scattering in the radiative transport equation and detector response.

14. An apparatus according to claim 13, arranged to work with the radiative transport equation combined with the detector response that has been simplified by the following assumptions:

- that each particle is scattered a limited number of times - that the scattering angle is limited

- that only particles that reach the detector in a direction substantially parallel to the detector normal will be detected.

15. An apparatus according to any one of the claims 12-13, wherein - the attenuation operator means (110, 120) is arranged to generate attenuation operator based on an attenuation map

- the scattering operator means (130) is arranged to generate the scattering operator based on the attenuation map, and

- the detector response operator means is arranged to model the detector response operator based on the properties of the detector.

16. An apparatus according to claim 15, further comprising attenuation map generating means (110) arranged to receive transmission data from the camera system (10) and to use these transmission data to generate the attenuation map.

17. An apparatus according to claim 15, further comprising an attenuation map generating means (110) arranged to generate the attenuation map on the basis on experience and/or prior knowledge.

18. An apparatus according to any one of the claims 15-17, wherein the attenuation operator means (110, 120) are arranged to receive said attenuation map from the attenuation map generating means (110) and to generate said attenuation operator from said attenuation map.

19. An apparatus according to claim 14, further comprising - means for making a virtual slicing of the object parallel to the detector plane - means for assigning a distance to the detector for the slice - means (160) for treating the detector response for the virtual slice as a 2D convolution with a convolution kernel dependent on the distance and direc- tion to the 2D detector - means for combining the treated detector responses a 3D convolution evaluated on the detector plane.

20. An apparatus according to claim 15, wherein the scattering operator means (130) is arranged to define the scattering operator as a 3D reconstruction with a kernel derived from the attenuation map.

21. An apparatus for high-resolution reconstruction of an object, arranged to use an iterative reconstruction method for the reconstruction, said apparatus comprising: - camera means (10) for receiving from at least one 2D detector information about particles emitted from the object, - an apparatus according to any one of the claims 14-20 arranged to receive camera parameters and transmission data from said camera means, and - model operator means for supplying said model operator to a device (90) ar- ranged to calculate a goodness of fit measure between the measured data and synthetic data generated by the model operator.

22. An apparatus according to claim 23, further comprising filtered backprojection means (100) for generating a first image on the basis of said transmission data re- ceived from the camera means (10)