Abstract:
A method for estimating a real function that describes a phenomenon occurring in a space of any dimensionality is disclosed. The function is estimated by taking a series of measurements of the phenomenon being described and using those measurements to construct an expansion that has a manageable number of terms. A reduction in the number of terms is achieved by using an approximation that is defined as an expansion on kernel functions, the kernel functions forming an inner product in Hilbert space. By finding the support vectors for the measurements one specifies the expansion functions. The number of terms in an estimation according to the present invention is generally much less than the number of observations of the real world phenomenon that is being estimated. In one embodiment, the function estimation method may be used to reconstruct a radiation density image using Positron Emission Tomography (PET) scan measurements.

Description:
This is a continuation-in-part application of U.S. patent application Ser. No. 08/726,452 filed Oct. 4, 1996, now U.S. Pat. No. 5,950,146. 
    
    
     The U.S. Government has a paid-up license in this invention and the right in limited circumstances to require the patent owner to license others on reasonable terms as provided for by the terms of Grant No. PHY 9512729 awarded by the N.S.F. 
    
    
     FIELD OF INVENTION 
     This invention relates to the fields of function approximation, regression estimation and signal processing. More particularly, the invention relates to practical methods for estimating functions of any dimensionality. This invention has particular utility in the reconstruction of images associated with Positron Emission Tomography (PET) scan data. 
     BACKGROUND OF THE INVENTION 
     Function estimation is a branch of mathematics having many real world applications. For example, function estimation may be applied to a finite number of noisy measurements (or “observations”) for the purpose of approximating a function that describes the phenomenon being measured. Existing techniques of function estimation generally use the expansion:                f        (   x   )       =       ∑     i   =   1     N            α   i            φ   i          (   x   )                   (   1   )                                
     where f(x) is the approximating function, φ i (s) are a plurality of a priori chosen functions, α i  are scalars, and “N” is the number of terms in the expansion. Accordingly, to estimate a function using equation (1), one must estimate the scalars α i . The equation may be readily applied to one or two dimensional data, however its application to data of three or more dimensions is limited by the “curse of dimensionality”. 
     The curse of dimensionality is well known among mathematicians. It refers to the exponential increase in complexity of the estimating function that occurs as the dimensionality of the data being estimated increases. For example, if a one dimensional curve may be estimated with reasonable accuracy by using 20 terms (N=20), then the estimation of a two dimensional surface with reasonable accuracy requires 20 2  terms (N=20 2 ). A three dimensional surface requires 20 3  terms, and so on. Thus, as the dimension of the function to be estimated increases, the number of terms required rapidly reaches the point where even the most powerful computers can not handle the estimation in a timely manner. 
     Positron Emission Tomography (PET) is a graphical method of recording the physiological processes of the body. A tracer labelled with a positron emitting nuclide is injected into the body, and an image of the distribution of radioactivity which follows is reconstructed from the measured positron emissions. 
     There are two classes of image reconstruction algorithm used in PET. Direct analytical methods because they require relatively short computation times. Images reconstructed by analytical methods tend to have low visual quality and display interference between regions of low and high tracer concentration. Iterative methods generally produce higher visual quality, but are computationally much more intensive. At present, the computation time of reconstruction using iterative methods is unacceptably long. 
     At present, analytical and iterative methods require the problem to be discretized at some stage. This means that the activity concentration is calculated at a finite number of points from a finite number of projections. The area to be reconstructed is subdivided into a mesh of pixels, and the activity concentration within each pixel is assumed to be uniform. The lines of response are rays of a finite width which depends on the detectors&#39; resolution. The “coincidence counts” for each line of response are proportional to the line integral of the tracer concentration. The Radon transform operator, R[f(x,y)] defines the integral of f(x,y) along a line, and analytical methods approximate the inverse Radon operator to determine the function f. 
     Filtered Back Projection (FBP) is an analytical PET image reconstruction method in which the emission data is first frequency filtered and is then back-projected into the image space. Back-projection of unfiltered emission data, p, is known as simple back-projection. Back-projection of p assigns a value which is proportional to the number of coincidence counts for a line of projection to each pixel through which that line passes. The resulting image is a star like pattern, with pixel values peaking at the position of the source. This means that nearby pixels will have non-zero values although no tracer is present at that location. The star-like image may be significantly reduced by initial frequency filtering of the emission data. The peak is principally composed of high-frequency components, while areas outside the peak are composed of low-frequency components. High-pass frequency filtering the emission data amplifies the high-frequency components relative to the low-frequency components, and the back-projection of the filtered projection data is of higher visual quality. 
     Iterative methods of reconstruction require that the emission data is discretized before solution of the reconstruction problem. A finite-dimensional algebraic system of linear equations are solved by some iterative method. Iterative methods of reconstruction are more flexible than analytical methods because weights or penalties which reflect the nature of the problem may be introduced into the set of linear equations. For example, additional constraints and penalty functions can be included which ensure that the solution has certain desired properties, such as requiring that values of neighboring pixels should be close to one another. 
     Iterative PET image reconstruction methods may be divided into three classes: Norm Minimizations, Maximum Likelihood Approaches, and Maximum A Posteriori Approaches. The prior art includes many Norm Minimization methods, although those which involve factorization of the matrix are usually unsuitable because of the size of the computational power required for solution. One set of norm minimization method, Algebraic Reconstruction Techniques, or Row Action Methods, are relaxation methods which require access to only one row of matrix per iteration. Maximum Likelihood Approaches (MLAs) maximize the conditional probability of observing the coincidence counts over all possible image vectors. MLA models which assumes the emission process has a Poisson distribution are considered to be the most appropriate to PET reconstruction, and are usually solved by the expectation maximization (EM) iterative method. Maximum A Posteriori Approach (MAPA) models maximize the conditional probability of the image vector given the measurement vector or observed coincidences. MAPA models require that the a priori probability distribution of the images is known, which has been considered a drawback. A recent comparison of the various iterative methods of PET image reconstruction concluded that a Maximum Likelihood Approach using EM methods best characterized the emission process and were the most promising direction for further PET image reconstruction research. 
     SUMMARY OF THE INVENTION 
     The present invention provides a method for approximating functions of any dimensionality by an estimating function having a manageable number of terms. The reduction in the number of terms is achieved by using the expansion:                f        (   x   )       =         ∑     i   =   1     l            a   i          K        (     x   ,     x   i       )           +   b             (   2   )                                
     where K(x,x i ) is a positively defined kernel function having two vector variables that define an inner product in some Hilbert space. The a i  are scalars, and the x i  are subsets of the measured data called support vectors (SVs). The number of terms required when approximating a function with this expansion depends on the complexity of the function being estimated and on the required accuracy of the approximation, but not on the dimensionality of the function being estimated. The invention allows the PET image reconstruction problem to be solved efficiently using a novel Norm Minimization approach. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a flow chart depicting the steps taken in using the support vector method for estimating a function according to the present invention. 
     FIG. 2 shows the geometry of a two dimensional Positron Emission Tomography (PET) scan machine. 
     FIG. 3 shows the basic geometry of the PET scan image reconstruction problem. 
     FIG. 4 is a block diagram of a system for estimating real functions according to the present invention. 
    
    
     DETAILED DESCRIPTION 
     The SV Method for Estimation of Real Functions 
     FIG. 1 is a flow chart depicting the steps taken in estimating a real function according to the present invention. The real function to be estimated will generally be descriptive of some physical phenomenon that occurs in an “n” dimensional space. Accordingly, the first step in estimating a real function is to take measurements of the phenomenon that is to be described by the estimated function (step  10 ). For purposes of illustration, the number of measurements is assigned the symbol l. Each of the l measurements is described by a vector “x” and a value “y”. The vector x i  describes the position within the n dimensional space at which the i th  measurement occurs, and therefore x is an n dimensional vector. The value y i  describes the value of the i th  measurement. Thus, the measured data is denoted as (x i ,y i ), (x 2 ,y 2 ) . . . (x i ,y i ). 
     Once the measurements have been acquired, the following functional is formed (step  12 ) and minimized (step  14 ):                  R        (   ω   )       =       C          ∑     i   =   1     l                   y   i     -     f        (       χ   i     ,   ω     )              l         +     (     ω   ,   ω     )         ,           (   3   )                        y   -     f        (     χ   ,   ω     )              l     =     {         0           if              |       y   -     f   (     x   ,   ω            &lt;   ɛ       ,                      y   -     f        (     x   ,   ω     )              -   ɛ         otherwise                   (   4   )                                
     (ω,ω) is the inner product of two vectors, C is some constant, and f(x,ω) is a set of functions that can be represented by:                f        (     x   ,   ω     )       =       ∑     i   =   1     ∞            ω   i            φ   i          (   x   )                   (   5   )                                
     The ε referred to in equations (3) and (4) is an error tolerance which is chosen a priori. It describes the maximum amount that the estimated function may deviate from the actual function described by the measured data. 
     The function that minimizes the functional described in (3) is: where α i *,α i ≧0 with α i α i =0; and (Φ(x i ),Φ(x)) is the inner product of two elements of Hilbert space, where Φ(x) is defined as Φ(x)=(Φ 1 (x),Φ 2 (x) . . . Φ k (x), . . . ).                f        (     x   ,   α   ,     α   *       )       =         ∑     i   =   1     l            (       α   i   *     -     α   i       )          (       Φ        (     x   i     )       ,     Φ        (   x   )         )         +   b             (   6   )                                
     Thus, in order to minimize the functional, coefficients α i * and α i  of equations (3) must be determined (step  16 ). This is done by solving the following quadratic optimization problem:                  W        (       α   *     ,   α     )       =         -   ɛ            ∑     i   =   1     l          (       α   i   *     +     α   i       )         +       ∑     i   =   1     l          y        (       α   i   *     -     α   i       )         -       1   2            ∑       i                 j     =   1     l            (       α   i   *     -     α   i       )          (       α   j   *     -     α   j       )          (       Φ        (     x   i     )       ,     Φ        (     x   j     )         )               ,           (   7   )                                
     subject to the constraints:                    ∑     i   =   1     l          (       α   i   *     -     α   i       )       =   0     ,     0   ≤     α   i       ,       α   i   *     ≤   C     ,     i   =   1     ,   …              ,   l           (   8   )                                
     With respect to the solution of equation (6), it is important to note that only some of the coefficients α i *-α i  differ from zero. Also, in equation (6) the vectors x i  are called Support Vectors (SV s ) and equation (6) is said to describe an expansion on these SVs. 
     In both equations (6) and (7), one can evaluate the inner products (Φ(x i ), Φ(x)) by using the general form of the inner product in Hilbert space. According to Hilbert space theory, to guarantee that a symmetric function K(u, v) has an expansion:                K        (     u   ,   v     )       =       ∑     k   =   1     ∞            a   k            φ   k          (   u   )              φ   k          (   v   )                   (   9   )                                
     with positive coefficients α k &gt;0, i.e. to guarantee that K(u, v) is an inner product in some feature space Φ, it is necessary and sufficient that the conditions 
      ∫K(u,v)g(u)g(v)dudv&gt;0  (10) 
     be valid for any non-zero function g on the Hilbert space (Mercer theorem). Therefore, equation (6) can be replaced with                             f        (     x   ,   α   ,     α   *       )       =         ∑     i   =   1     l            (       α   i   *     -     α   i       )          K        (     x   ,     x   i       )           +   b               (   11   )                                
     where the inner product (Φ(x i ), Φ(x)) is defined through a kernel K(x i ,x). In which case, to find coefficients α i * and α i  one has to maximize the function:                W        (       α   *     ,   α     )       =         -   ɛ            ∑     i   =   1     l          (       α   i   *     +     α   i       )         +       ∑     i   =   1     l          y        (       α   i   *     -     α   i       )         -       1   2            ∑       i                 j     =   1     l            (       α   i   *     -     α   i       )          (       α   j   *     -     α   j       )          K        (       x   i          x   j       )                       (   12   )                                
     subject to the constraints of equation (8). 
     Coefficients α i * and α i  are scalars, and therefore equation (11) may be rewritten as:                f        (   x   )       =         ∑     i   =   1     l            a   i          K        (     x   ,     x   i       )           +   b             (   13   )                                
     where K(x,x i ) is a positively defined kernel function having two vector variables that define an inner product in some Hilbert space. The a i  in equation (13) are scalars, and the x i  are subsets of the measured data called support vectors. The number of terms required when approximating a function with the expansion of equation (13) depends on the complexity of the function being estimated and on the required accuracy of the approximation, but not on the dimensionality of the function being estimated. 
     Use of the Support Vector Method to Solve Linear Operator Equations 
     The SV method may be used for solving linear operator equations 
     
       
         Af(t)=F(x)   (14) 
       
     
     where operator A realizes a one-to-one mapping from a Hilbert space E 1  into a Hilbert space E 2 . 
     Instead of the function F(x) on the right-hand side of (14) we are given measurements of this function, generally with errors. 
     
       
         (x 1 ,F 1 ), . . . ,(x l ,F l ).  (15) 
       
     
     It is necessary to estimate the solution of equation (14) from the data (15). 
     If the solution to (14) belongs to the set of function f(t, w)∈Φ linear in some flattening space {φ 0 (t), . . . ,φ N (t), . . . }:                  f        (     t   ,   w     )       =         ∑     r   =   0     ∞            w   r            φ   r          (   t   )           =     (     W   ,     Φ        (   t   )         )         ,           (   16   )                                
     where we denote 
     
       
         W=(w 1 , . . . , w N , . . . ), 
       
     
     
       
         Φ(t)=(φ 0 (t), . . . ,φ N (t), . . . ). 
       
     
     To use the SV technique we solve this equation by minimizing the functional                D        (   f   )       =       C          ∑     i   =   1     l                   A                 f     -     F   i            ɛ         +     γ        (     W   ,   W     )                 (   17   )                                
     where the regularization term is defined by the same operator than defines the equation. 
     The operator A maps this set of functions into the set of functions                  F        (     x   ,   w     )       =       A                   f        (     x   ,   w     )         =         ∑     r   =   0     ∞            w   r        A                     φ   r          (   t   )           =         ∑     r   =   0     ∞            w   r            ψ   r          (   x   )           =     (     W   ,     Ψ        (   x   )         )             ,           (   18   )                                
     which are linear in another flattening space 
     
       
         Ψ(x)=(ψ 0 (x), . . . ,ψ N (x), . . . )  (19) 
       
     
     
       
         where  
       
     
     
       
         ψ r (x)=Aφ r (t).  (20) 
       
     
     To minimize functional (17) in the set of functions f(t,w), it is sufficient to find the vector coefficients W, which minimizes the functional                  D   1          (   F   )       =       C          ∑     r   =   0     l                   F        (     x   i     )       -     F   i            ɛ         +     (     W   ,   W     )               (   21   )                                
     in an image space that is in the space of functions F(x i ,w). The generating kernel in the image space is defined as                K        (       x   i     ,     x   j       )       =       ∑     r   =   0     ∞              ψ   r          (     x   i     )              ψ   r          (     x   j     )                   (   22   )                                
     and the cross-kernel function is defined as                κ        (       x   1     ,   t     )       =       ∑     r   =   0     ∞              ψ   r          (     x   i     )                φ   r          (   t   )       .                 (   23   )                                
     The problem of finding the solution to the operator equation (finding the corresponding vector of coefficients W) is equivalent to the problem of finding vector W for linear regression function (18) in the image space using measurements (15). 
     This regression problem may be solved using the quadratic optimization SV technique, maximizing the objective function (12) subject to constraint (8). The kernel (22) is used to find both the support vectors x i , i=1, . . . , N and the corresponding coefficients α i *-α i  that define the vector W for SV regression approximation              W   =       ∑     i   =   1     N            (       α   i   *     -     α   i       )          ψ        (     x   i     )                   (   24   )                                
     It is sufficient to use the standard quadrative optimization SV technique. The coefficients W with the regression in image space define the approximation to the desired solution, so inserting these coefficients into (16) yields                f        (     t   ,   α   ,     α   *       )       =       ∑     i   =   1     N            (       α   i   *     -     α   i       )            κ        (       x   i     ,   t     )       .                 (   25   )                                
     Solving of the operator equation thus involves using the cross-kernel function as an expansion on the support vector. 
     In order to solve the linear operator equation using the SV method, the following steps are taken: 
     1. Define corresponding regression problem in image space. 
     2. Construct the kernel function K(x i ,x j ) for solving the regression problem using the SV method. 
     3. Construct the cross-kernel function k(x i ,t). 
     4. Use the kernel function K(x i ,x j ) to solve the regression problem by the SV method (i.e. to find the support vectors x 1 *, i=1, . . . , N and the corresponding coefficients β i =(α i *−α i , i=1, . . . N). 
     5. Use these supports vectors and corresponding coefficients to define the solution: In these five steps the first three steps (constructing the regression problem, constructing the kernel in the image space, and constructing the corresponding cross-kernel function) reflect                f        (   t   )       =       ∑     r   =   1     N            β   r          κ        (       x   r     ,   t     )                   (   26   )                                
     the singularity of the problem at hand (they depend on the operator A). The last two steps (solving the regression problem by the SV method and constructing the solution to the desired problem) are routine. However, a given operator equation must obtain explicit expressions for both the kernel function in image space and the corresponding cross-kernel function. 
     The solution of ill-posed operator equations requires additional steps. When solving the operator equation where only the approximation F l  of the right-hand side 
     
       
         Af=F  (27) 
       
     
     of (27) is given, the functional must be minimized in a set of functions {f}. In this 
     
       
         W(f)=||Af−F t || 2 +γ l Ω(f) 
       
     
     functional the term Ω(f) is the regularization functional and the parameter γ l  is the regularization constant. The constant γ l  is chosen using the discrepancy principle. If, for sufficiently large l, the accuracy of approximating function F l  obtained from the data does not exceed ε, then regularization function Ω(f) must be minimized subject to 
     
       
         ||Af−F l ||≦ε  (28) 
       
     
     By using a loss function insensitive to ε, the SV method of solving operator equation applies this principle. For sufficiently large C, the SV method minimizes the regularization functional, which is the norm of the vector of coefficients of the linear function in flattening space subject to constraint (28). 
     Where detailed information on accuracy of the approximation is available, the regularization process may be better controlled. Data describing the right-hand side of (27) may be obtained simultaneously with an estimate of the accuracy of obtaining data at any specific set of points. In other words, the operator equation is solved given triples 
     
       
         (x 1 , F 1 , ε 1 ), . . . ,(x l , F l , ε l ).  (29)  
       
     
     Thus, using various ε-insensitivities for various points (vectors) x, the regularization process may be better controlled. 
     Use of the Support Vector Method to Estimate Functions Generated by Positron Emission Tomography 
     Positron Emission Tomography (PET) is medical diagnostic technique which involved the reconstruction of radio activity within the body following the injection or inhalation of a tracer labeled with a positron emitting nuclide. During the disintegration of the radioactive nuclei collected in the body, positrons are emitted. These positrons collide with nearby electrons resulting in the annihilation of the electron and positron and the emission of two gamma rays in opposite directions. FIG. 2 shows the detection of the gamma rays on opposite sides of the radioactive source. From each point within the source  200 , the gamma ray pair can be emitted in any direction. Two-dimensional PET, however, measures only rays that belong to a fixed plane  202 . In this plane if gamma ray  204  hits a first detector  206  and then within a small time interval another gamma ray  208  hits the opposite dectector  210 , then it is known that an emission must have originated from a point somewhere along the line joining these two detectors, the line of response  212 . This event is known as a coincidence. The total number of coincidences for this pair of detectors is proportional to the integral of the tracer concentration along line of response  212 . In order to obtain information about the tracer distribution in the plane, PET devices include a large number of detector pairs with lines of response at many different angles. The set of all detector pairs whose lines of response are at a given angle μ form a μ-projection. A sinogram contains the set of all projections. At present, a typical PET device includes between 100 and 300 projection angles, μ j  with each projection having between 100 and 200 lines of response m i . This gives between 10,000 and 60,000 lines of response, each with corresponding recorded number of coincidences p(m k ,μ k ). Thus, the data includes l triplets m k ,μ k ,p(m k ,μ k ) called observations in each sinogram. The object is to reconstruct the density of nuclear concentration within a given plane of the body using the sinogram. 
     FIG. 3 shows a line or response  300  inclined at angle μ to the y axis and at distance m from the origin  302 . The detector circle  304  is assumed to have a radius=1. Coincidence count p(m,θ) is proportional to integral of concentration function f(x,y) along a line  306  defined by a pair m,μ. The Radon operator R defines the integral of f(x,y) along any line                  R        [     f        (     x   ,   y     )       ]       =         ∫     -       1   -     m   2               1   -     m   2                  f        (         m                 cos                 μ     +     u                 sin                 μ       ,       m                 sin                 μ     -     u                 cos                 μ         )               u         =     p        (     m   ,   μ     )           ,           (   30   )                                
     where coordinates x, y along the line are defined by equations 
     
       
         x=m cos μ+u sin μ 
       
     
     
       
         y=m sin μ−u cos μ  (31) 
       
     
     and the position of the line is defined by parameters 
     
       
         −1&lt;m&lt;1, 0≦μ≦π.  (32) 
       
     
     The interval of integration is defined by 
     
       
         −a m   =−{square root over (1+L −m 2 +L )}≦u≦{square root over ( 1+L −m 2 +L )}=+a m   (33) 
       
     
     The principal advantage of the use of the Radon transform operator to find the function f(x,y) that satisfies (30) given the function p(m,θ) is that under wide conditions inverse operator 
     
       
         R −1 R[f(x,y)]=R −1 p(m,μ)=f(x,y)  (34) 
       
     
     exists. In order words, the solution to the Radon equation exists. However, the finding of this solution is an ill-posed problem. 
     The object is to find the solution to this ill-posed problem in a situation where function p(m,θ) is defined by its values p k  in a finite number l of points m k ,μ k , k=1, . . . l. Moreover, the data points are corrupted by some random noise 
     
       
         p k =p(m k ,μ k )  (35) 
       
     
     The object is to solve the PET Radon equation (30) given measurements 
     
       
         (p 1 , m 1 , μ 1 ), . . . , (p l , m l , η l )  (36) 
       
     
     which is a stochastic ill-posed problem. 
     In order to find the solution to the operator equation 
     
       
         Af(t)=F(x)   (37) 
       
     
     using approximation F l (x) instead of exact right hand side F(x) of (37), the regularized functional 
     
       
         R=||Af(t)−F * (x)|| 2 +γΩ(f),   (38) 
       
     
     must be minimized in a given set of functions {f(t)}, where γ&gt;0 is some regularization constant and Ω(f) is a regularizing functional. 
     In the PET problem, where a finite number of measurements are given, the following functional may be constructed: The set of functions in which one looking for the solution may be selected from piece-wise constant or piece-wise linear functions {f(x,y)}.                R        (   f   )       =         ∑   k   l            (       P   k     -       ∫     -   1     1            f        (           m   k        cos                   μ   k       +     u                 sin                   μ   k         ,         m   k                   sin                   μ   k       -     u                 cos                   μ   k           )               u           )     2       +     γΩ        (   f   )                 (   39   )                                
     The residual method may be applied to solve an integral equation with the right-hand side of the equation approximately defined. If the linear operator equation (31) is to be solved with approximation F l (x) instead of exact equation F(x), and information about the accuracy of the approximation 
     
       
         Δ=||F(x)−F l (x))||.  (40)  
       
     
     is available, the residual method suggests the choice of f γ (t) which minimizes the functional ||f|| 2  and satisfies the constraint 
     
       
         ||Af γ (t)−F l (x))||≦Δ  (41) 
       
     
     For PET image reconstruction, the exact value of Δ cannot be evaluated, because it is the result of measurement of a random event in which stochastical accuracy is characterized by the variance. The random value of the number of coincidences along any line of response may be characterized as 
     
       
         ε k ={square root over (p(m k +L ,μ k +L ).)}  (42) 
       
     
     For PET image reconstruction the functional ||f|| 2  is minimized subject to constraints                       p   k     -       ∫     -     a   m         +     a   m                  f   γ          (           m   k        cos                   μ   k       +     u                 sin                   μ   k         ,         m   k        sin                   μ   k       -     u                 cos                   μ   k           )               u                ≤     δɛ   k             (   43   )                                
     where δ&gt;0 is some constant. As described above, the SV method with a ε-insensitive loss function (with different ε i  for different vectors x i ) is an appropriate method for solving such problems. 
     Prior art solution methods of solving the PET image reconstruction problem find the solution in a set of piece-wise constant functions. For this purpose the n×n=N pixel space was introduced where in any pixel the value of function was considered to be constant. If l is the number of measurements, the integral equation (30) was approximated by the algebraic equation 
     
       
         Ax=b, x≧0,   (44) 
       
     
     where A∈R l×N  was a known matrix, x∈R N  was a vector defined the values of the approximating function in the set of pixels, and b∈R l  was a vector that defined the number of coincidences along lines of response. 
     The regularized method for solution of equation (44) proceeded by minimizing the functional 
     
       
           R= ( b−Ax ) ( b−Ax ) T +Ω( x ),  (45) 
       
     
     where the regularizer Ω(x)=(x,x) was used. Other regularizers were also used. 
     Use of the strong residual principle for solving this problem led those in the prior art to choose a regularization parameter γ for which the solution x* of the optimization problem (44) minimized the functional (x,x) and satisfied the constraints                         b   k     -       ∑     i   =   1     l            a     i   ,   k            x   k   *                ≤     ɛ   k       ,           (   46   )                                
     where b x  is coordinate of vector b, b i * is coordinate of vector x*, and a i,k  is element of matrix A. 
     The main problem with this prior art method of solving the PET image reconstruction problem is the computational power required for solving an equation this large. As described above, M may be as large as 60,000 observations and the number of parameters to be estimated N may also be approximately 60,000 (N=256×256). The prior art methods of solving the PET image reconstruction problem considered various methods of solving this equation, but to date only the SV method does not require the reduction of the problem to solving of the system of linear algebraic equations with a very large number of variables. 
     In accordance with the present invention, solutions to the PET image reconstruction problem may be sought in the set of two dimensional tensor-product spline functions built from one dimensional splines with knot sequence            {     t   i     }       t   =   i     N     .                          
     The knot spacing may be chosen as small as the resolution of the application requires; increasing the resolution in this way has very little effect on the number of support vectors, and hence on the difficulty of obtaining the solution. To simplify the formulas, piece-wise constant functions are considered. The construction of higher order splines is analogous. The desired function may be approximated by the expression                f        (     x   ,   y     )       =       ∑     i   ,   j              ψ     i   ,   j            θ        (     x   -     t   i       )            θ        (     y   -     t   j       )                   (   47   )                                
     where 
     ψ i,j  are the expression coefficients and θ is the Heaviside step function. 
     Transforming this function with the Radon operator R yields the corresponding regression problem in image space                p        (     m   ,   u     )       =       ∑     i   ,   j              Ψ     i   ,   j            φ        (     m   ,   μ   ,     t   i     ,     t   j       )                   (   48   )                                
     where we have set 
     
       
         p(m,μ)=(Rf)(m,μ) 
       
     
     
       
         φ(m,μ)=(Rθ(−t i )θ(−t j ))(m,μ).  (49) 
       
     
     The problem of solving the PET image reconstruction equation in the set of piece-wise constant functions is reduced to the problem of regression approximation in the image space using the data          {       p   i     ,     m   i     ,     μ   i       }            l     i   =   1       .                            
     The solution requires finding the constants ψ i,j  that satisfy the conditions                 p   k     =       ∑     i   ,   j              ψ     i   ,   j            φ        (       m   k     ,     μ   k     ,     t   i     ,     t   j       )                  ≤     δε   k                            
     and has a minimal norm. In the present invention, this is done by minimizing (21) to obtain the solution of the PET image reconstruction problem. In order to accomplish this, we must construct two functions: the kernel function in image space                K        (       m   i     ,     μ   i     ,     m   j     ,     μ   j       )       =       ∑     α   ,   β              φ        (       m   i     ,     μ   i     ,     t   α     ,     t   β       )            φ        (       m   j     ,     μ   j     ,     t   α     ,     t   β       )                   (   50   )                                
     and the cross-kernel function                κ        (       m   i     ,     μ   i     ,   x   ,   y     )       =       ∑     α   ,   β              φ        (       m   i     ,     μ   i     ,     t   α     ,     t   β       )            φ        (     x   -     t   α       )            φ        (     y   -     t   β       )                   (   51   )                                
     These expressions may easily be evaluated analytically. 
     To solve the PET image reconstruction problem using the SV technique, two steps remain. First, the kernel-function (50) is used to obtain the support vectors (m k ,μ k ), k=1, . . . , N and corresponding coefficients α* k −α k , k=1, . . . N. Second, the cross-kernel function (51) is used to obtain support vectors and coefficients that define the desired image reconstruction approximation                f        (     x   ,   y     )       =       ∑     k   =   1     T            (       α   k   *     -     α   k       )          κ        (       m   k     ,       μ   k     ;   x     ,   y     )                   (   52   )                                
     Note that for any support vector (m k , μ k ) in image space there is a corresponding line of response in pre-image space defined by expression 
     
       
         x=m k cos μ k +usin μ k   
       
     
     
       
         y=m k sin μ k −ucos μ k    
       
     
     
       
         {square root over (1−+L m k   2 +L )}≦u≦{square root over (1+L −m k   2 +L )}  (53) 
       
     
     Therefore, in pre-image space the expansion of the function using support vectors (m k , μ k ), k=1, . . . , N is equivalent to the expansion of the desired solution on the lines of response. 
     FIG. 4 is a block diagram of a system for estimating real functions according to the present invention. As mentioned above, the function being estimated will generally describe an “n” dimensional phenomenon of which “l” measurements have been taken. In the system of FIG. 5, these measurements are provided to a function estimator  400  through a data source  402 . The data source  402  may be a set of PET detectors. The estimator operates on the measured data to generate the estimating function. 
     The data source and function estimator may take many forms. For example, the data source may be the measuring instrument used to generate the measured data, or it may merely be a storage device, such as a disk drive, where measured data is stored. The function estimator may be computer, such as a personal computer, or it may be a single circuit board, such as an expansion board for a personal computer. 
     In any event, the function estimator includes an input device  404  that receives the measurements from the data source. An example of an input device is an RS-232 port (commonly available on many personal computers). The input device, in turn, relays the measurements to a processing unit  406 , where the estimation process of forming a functional, minimizing the functional, and determining the coefficients of the minimizing function is implemented. 
     In a preferred embodiment, the processing unit includes a memory  408 , one or more stored programs  410 , and a processor  412 . The memory is used to store measurement data received through the input device. It may be, for example, the Random Access Memory (RAM) of a computer. The stored programs are those programs necessary for implementing the estimation process. They may be stored in the Read Only Memory (ROM) of a computer. The processor, which may be the Central Processing Unit (CPU) of a computer, implements the estimation process according to instructions set forth in the stored programs. 
     In a preferred embodiment, the function estimator includes an output device  414 . The output device is used to provide information to the user of the function estimator by, for example, providing a plot of the estimated function superimposed on a plot of the measured data. An example of an output device suitable for use with the invention is a computer display monitor. 
     Several embodiments for the invention have been described in detail. However, it should be obvious that there are many further variations evident to persons skilled in the art. The invention is more particularly defined in the appended claims.