Patent Publication Number: US-2009221932-A1

Title: Microwave temperature image reconstruction

Description:
TECHNICAL FIELD OF THE INVENTION 
     The present invention relates to a microwave temperature image reconstruction system. More specifically, the present invention relates to the use of an algorithm which shows itself particularly adapted to temperature image reconstruction. 
     BACKGROUND OF MICROWAVE IMAGING 
     The large potential of temperature imaging has been recognized for long. For example in the medical community such data is known to potentially improve a variety of diagnostic and interventional procedures. In the case of industrial applications, thermal information has shown to predict and prevent failures of electrical circuits and equipment. One potentially interesting imaging modality for temperature imaging is microwave imaging. Due to several practical drawbacks, microwave imaging has hardly been able to enter the industrial or medical end-user market though. One of these drawbacks lies in the difficult image reconstruction task which contrasts with the rather simple data acquisition process of a potential microwave imaging system. 
     Non-Destructive Evaluation (NDE) has been one of the engineering disciplines which mostly revolutionized diagnostic techniques in industry and in medicine during the last decades. MR (Magnetic Resonance), CT (Computerized Tomography), US (Ultra-Sound), and other NDE devices are being standard tools for a wide range of diagnostic fields. Furthermore they are currently changing significantly medical surgery, as their capability of visualizing intra-body information enables surgeons to minimize the invasiveness of their interventions. Even though NDE techniques are by themselves expensive, they are potentially even interesting from a financial point of view, as they can significantly decrease the expensive hospitalization time of patients. Similar arguments apply to industrial NDE which was able to bring quality assurance and failure prediction to an impressive performance. 
     Several NDE modalities are currently being used, including MR, CT, and US imaging. Each of these modalities images a particular physical quantity, such as MR can image proton density, CT images x-ray absorption coefficients, or US visualizes reflection coefficients for ultrasound waves. Therefore each of these modalities has applications for which they are mostly adapted, while being excluded from other applications due to their incapability to image the appropriate physical quantity. As a result, the development of new NDE imaging modalities which are able to image additional physical quantities is a main branch of engineering and fundamental research. One of these quantities, which have been recognized to be particularly valuable for a wide range of medical and industrial applications, is intra-body temperature. Currently, just MR imaging is capable to perform such measurements non-invasively, though being prohibitively expensive for a wide range of promising applications. As an example, cancer development is known to be accompanied with a local temperature increase of 1-2° C. Nevertheless, the temperature resolution and mere cost of MR temperature imaging prohibits its use in the important field of breast cancer detection. 
     As an alternative to MR temperature imaging, passive microwave temperature monitoring is a promising avenue. Microwave temperature imaging devices are in comparison to MR devices significantly smaller, cheaper, and less complex to implement. Nevertheless, their current research implementations are too premature to be produced and marketed on large scales. One of the main drawbacks lies in the fact that the data collected by microwave imaging devices is very sparse (only around 7 independent measurements) in comparison to the temperature resolution and spatial resolution one would like to obtain for the final temperature profiles. This implies that the image reconstruction is particularly challenging for such devices. 
     SUMMARY OF THE INVENTION 
     The present invention proposes a solution to this reconstruction problem. It uses “large” over-complete dictionaries which, on the one hand, reflect the variability of potential temperature profiles and, on the other hand, incorporate as much a-priory information as possible in order to provide reliable image reconstruction. Together with state-of-the-art algorithms for data decomposition using over-complete dictionaries, temperature profiles can be reconstructed reliably. 
     The system according to the invention uses passive microwave imaging for temperature imaging, where “passive” refers to the fact that the imaging system does not emit any radiation into the investigated object. Rather, it senses and measures the natural thermal radiation emitted from this object, which contrasts to active and scattering-based microwave imaging devices. 
     The above and other objects, features and advantages of the present invention will become apparent from the following description when taken in conjunction with the accompanying drawings which illustrate preferred embodiments of the present invention by way of example. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  shows 3 B-splines of a set of 100 B-splines which cover regularly the imaging space and whose weighted sum constitutes the temperature profile reconstruction. 
         FIG. 2  shows three results of current state-of-the-art microwave temperature profile reconstruction by Tikhonov regularization and its limitations. The reference and reconstructed profiles are shown on the left side, while the resulting expansion coefficients are presented in the right column. 
         FIG. 3  shows the basic building blocks necessary for a microwave temperature monitoring device with a reconstruction procedure as the one disclosed herein. 
         FIG. 4  shows a schematic sketch of the functional building blocks necessary to realize the disclosed invention. 
         FIG. 5  shows the functional building blocks for data decomposition with over-complete dictionaries as used in the disclosed invention. 
         FIG. 6  presents a sketch of the building blocks for a PC-based implementation of the disclosed invention. 
         FIG. 7  shows a more detailed view of the PC-based sample implementation of a microwave temperature imaging device. 
         FIG. 8  shows cubic B-splines at scale 0 and scale 9 covering regularly the imaging space. These functions were part of the over-complete basis function set, used in the sample implementations of the disclosed invention. 
         FIG. 9  presents three reconstruction examples for the PC-based sample device. On the left side, the reconstructed and reference profiles are drawn. On the right side, the determined optimal expansion coefficients are shown. These coefficients correspond to the reconstruction profiles shown on the left. 
         FIG. 10  shows a sketch of the functional building blocks for the embedded sample implementation of the disclosed invention. 
         FIG. 11  shows a more detailed view of the embedded sample implementation of a microwave temperature imaging device. 
         FIG. 12  shows the flow chart of the FOCUSS algorithm as used for the temperature profile image reconstruction on the embedded sample implementation of the disclosed invention. 
         FIG. 13  shows three reconstruction results for the embedded sample device. The reconstructed and reference profiles are drawn on the left side, while the optimal expansion coefficients are drawn on the right side. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Passive microwave image reconstruction aims to solve the following equation for the unknown temperature profile, noted T(x): 
     
       
         
           
             
               
                 
                   
                     
                       b 
                       i 
                     
                     = 
                     
                       
                         γ 
                         i 
                       
                       · 
                       
                         
                           ∫ 
                           Ω 
                         
                          
                         
                           
                             
                               
                                 W 
                                 i 
                               
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                             · 
                             
                               T 
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                           
                            
                           
                               
                           
                            
                           
                              
                             x 
                           
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                    
                   
                     
                       with 
                        
                       
                         
 
                       
                        
                       i 
                     
                     = 
                     1 
                   
                   , 
                   2 
                   , 
                   … 
                    
                   
                       
                   
                   , 
                   M 
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   1 
                 
               
             
           
         
       
     
     where i indexes the M measurements of the system at frequencies f i . b i  are the passively measured brightness temperatures sensed by the multi-frequency imaging antenna at different frequencies f i , as described in [2, 9]. γ i  are the antenna efficiencies at frequencies f i , W i (x) are the spatial weighting functions at frequencies f i , and Ω is the sensing space of the microwave antenna. The spatial variable x may be a 1D, 2D, or 3D spatial vector, depending on the dimension of the reconstruction problem. In the 1D case (just one multi-frequency antenna), the number of measurements M lies at around 5 to 7, i.e. the antenna senses the brightness temperatures b i  at 5 to 7 different frequencies f i . Therefore we have just 5 to 7 measurements to recover the whole 1D temperature profile T(x). 
     The weighting function W i (x) weights the brightness temperature contributions from different spatial positions x. It therefore accounts for the absorption and scattering process of the temperature dependent thermal radiation from any point inside the investigated body to the passive measurement antenna. 
     Algorithms currently being used for the microwave image reconstruction of eq. 1 are presented in [1, 2]. Both approaches attempt to counter-balance the problem of having just few radiometric measurements available for the reconstruction, by giving a very restrictive prior on the temperature profile shape one would like to recover. In [1] and in the case of a 1D profile, the profile is assumed to be a single Gaussian: 
     
       
         
           
             
               
                 
                   
                     
                       T 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         T 
                         0 
                       
                       · 
                       
                         exp 
                          
                         
                           ( 
                           
                             - 
                             
                               ( 
                               
                                 
                                   x 
                                   - 
                                   
                                     x 
                                     0 
                                   
                                 
                                 σ 
                               
                               ) 
                             
                           
                           ) 
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   2 
                 
               
             
           
         
       
     
     with T 0  being the amplitude of the Gaussian, x 0  its position and σ its variance. The reconstruction algorithm determines the parameters, T 0 , x 0  and σ, which explain optimally the measurements b i  through eq. 1. The problem with this approach lies in the very strong Gaussian assumption of eq. 2. This assumption is too restrictive to deal with a large number of real world temperature profiles. 
     A more flexible approach to temperature image reconstruction has been proposed in [2]. There, a small sum of Chebyshev polynomials is assumed to build the temperature profile T(x): 
     
       
         
           
             
               
                 
                   
                     
                       T 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                        
                       
                           
                       
                        
                       
                         
                           t 
                           k 
                         
                         · 
                         
                           
                             h 
                             k 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   3 
                 
               
             
           
         
       
     
     where the h k (x) are the Chebyshev polynomials, and the t k  are the respective weighting coefficients. K is the number of polynomials in the sum which is assumed to result in a good approximation of the real temperature profile T(x) by eq. 3. Using eq. 3, eq. 1 can be re-written in matrix form: 
         b=A·t,   Eq. 4 
     with 
     
       
         
           
             
               A 
               
                 i 
                 , 
                 k 
               
             
             = 
             
               
                 γ 
                 i 
               
               · 
               
                 
                   ∫ 
                   Ω 
                   
                       
                   
                 
                  
                 
                   
                     
                       
                         W 
                         i 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     · 
                     
                       
                         h 
                         k 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                   
                    
                   
                       
                   
                    
                   
                     
                        
                       x 
                     
                     . 
                   
                 
               
             
           
         
       
     
     This notation leads quite naturally to matrix inversion and Tikhonov regularization as used in [2] to solve the microwave reconstruction problem with Chebychev polynomials or potentially other functions. This approach though has still some important shortcomings which will be illuminated in the following paragraphs, and which motivated the disclosed invention. 
     Let us shortly assume an example where the Chebychev polynomials of [2] are replaced by 100 B-splines of order 3 (cubic B-splines) at one single scale and covering the whole imaging depth. The sum of eq. 3 with 100 B-splines would be able to represent a very large variety of temperature profiles T(x). In  FIG. 1 , three of the 100 B-splines are drawn. Using all these 100 basis functions to reconstruct from the, let&#39;s say, 7 microwave measurements a temperature profile is not trivial, as one would try to solve a system of 7 equations for 100 unknowns t k  (eq. 4 with A being a 7×100 matrix, i.e. an ill-posed problem). Such a system is heavily underdetermined, which means that a large or probably even infinite number of solutions exists. In [2], Tikhonov regularization is proposed in order to deal with the ill-posed character of the resulting reconstruction problem. It applies a smoothness constrained by imposing a minimal l 2 -norm on the coefficients t, i.e. solving the following minimization problem in the least-square sense for t: 
       min∥A·t−b∥+λ·t·t T ,  Eq. 5 
     where λ weights the importance of the smoothing constraint when solving eq. 1 for the optimal temperature profile, T(x). On the contrary to eq. 4, eq. 5 is not an ill-posed problem, as long as λ stays sufficiently big. The least-square solution to eq. 4 is written: 
         t =( A   T   ·A+λ·I ) −1   ·A   T   ·b.   Eq. 6 
     Here, I is the identity matrix. One can see that when λ gets too small and the matrix A T ·A is ill-conditioned (which is in general the case for microwave temperature imaging) the inverse in eq. 6 does not exist. 
     In  FIG. 2 , we show some reconstruction results using the algorithm of [2] with 100 B-splines at one scale covering regularly the imaging space and which do not correspond to an over-complete set of basis functions. The reconstruction results are represented in broken lines and the references in full lines. One sees easily that this approach does not always lead to satisfactory results. The regularization term over-smoothes the reconstruction result when the temperature profile to be recovered is not smooth enough itself. On the other, it&#39;s not possible to decrease the smoothing parameter, λ, beyond some minimum, as this would result again in inverting an ill-conditioned matrix A T ·A+λ·I. On the right hand side of  FIG. 2 , one can also observe that this state-of-the-art approach to profile reconstruction results in a very large weighted sum of B-spline functions, as almost all K coefficients t k  are significantly larger than zero. In other words, Thikhonov regularization is neither appropriate for the described microwave image reconstruction problem. These arguments describe the motivations for the disclosed invention. They can be summarized as follows:
         1) The reconstruction algorithm has to deal with the ill-posed character of the reconstruction problem, while being applicable to a wide range of possible real-world temperature profiles.   2) The algorithm must not be dominated by any regularization term provoking over-smoothing or other artifacts to appear.       

     The general solution to these problems as disclosed in this invention can be summarized as follows:
         1) A large number K of basis functions, building an over-complete dictionary, has to be used in order to represent the variety of possible real-world temperature profiles with a small number N of these basis functions. As an example, the basis functions might be B-splines at different scales and positions. B-splines at several scales result in an over-complete set of basis functions. The fact that a small sum of N basis functions (small in the sense of N&lt;&lt;K) is sufficient for the reconstruction avoids over-smoothing of the resulting temperature profile estimation.   2) A regularized reconstruction algorithm has to be used which selects the optimal set of N basis functions out of the over-complete set of K initial basis functions (the over-complete dictionary) and which determines their weights to reconstruct the temperature profile. The regularization has to be as not to over-smooth the reconstruction result as shown to happen with e.g. Tikhonov regularization.       

     In conclusion, the disclosed invention describes a temperature image reconstruction procedure which, on the one hand, is capable to deal with the sparseness of the acquired radiometric data (just around M=7 measurements), and, on the other hand, succeeds in reconstructing a wide range of real-world temperature profiles for real-world medical and industrial applications without artifacts from the regularization term. 
     As will be shown in more detail in the following paragraphs, the disclosed invention maps the reconstruction problem into a well known formulation of decomposing signals with over-complete dictionaries and uses standard algorithms, such as matching pursuit [3], orthogonal matching pursuit [4], basis pursuit [5], high-resolution matching pursuit [6], etc. to solve the resulting equivalent problem. 
     The disclosed invention applies equivalently to 1D (profiles), 2D (images), or 3D (volumes) temperature image reconstruction. 
     Detailed Description of the Preferred Embodiments 
     The disclosed invention for radiometric temperature image reconstruction has to be implemented on a completely functional microwave temperature imaging device. Several design implementations of such a system are possible, all of which share common characteristics though. In  FIG. 3 , the functional building blocks common to all these systems are presented. First, we need one or several passive multi-frequency microwave antennas  1  for sensing the brightness temperatures, next a radiometer  2  has to interpret the signal and translate the sensed microwave radiation power into units of voltage. An analogue-to-digital converter  3  then digitizes the analogue voltage, so that the reconstruction unit  4 , basically any digital signal processing device, can reconstruct the intra-body temperature profile which caused the sensed microwave radiation. As a final step, the temperature profiles are visualized for inspection by the medical doctor or industrial investigator on a visualization unit, such as a computer screen  5 . The disclosed invention is physically situated in the reconstruction unit  4  and describes how temperature profiles can reliably be reconstructed from the sparse microwave radiation measurements. This reconstruction procedure will be described next. 
     Microwave image reconstruction aims to solve the following integral equation: 
     
       
         
           
             
               
                 
                   
                     
                       b 
                        
                       
                         ( 
                         f 
                         ) 
                       
                     
                     = 
                     
                       
                         γ 
                         f 
                       
                       · 
                       
                         
                           ∫ 
                           Ω 
                           
                               
                           
                         
                          
                         
                           
                             
                               W 
                                
                               
                                 ( 
                                 
                                   f 
                                   , 
                                   x 
                                 
                                 ) 
                               
                             
                             · 
                             
                               T 
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                           
                            
                           
                              
                             x 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   7 
                 
               
             
           
         
       
     
     where b(f) is the measured brightness temperature sensed by the passive imaging antenna at frequency f, γ f  is the antenna efficiency at frequency f, W(f,x) is the spatial weighting function at frequency f, Ω is the sensing space of the microwave antenna, and T(x) is the spatial temperature distribution the reconstruction algorithms aim to recover. The weighting functions W(f,x) might be obtained from a analytical model, as presented in [2] for the 1D case, or by electromagnetic simulations. The spatial variable x may be a 1D, 2D, or 3D vector, depending on the dimension of the reconstruction problem. Eq. 1 is known as a Fredholm equation of the first kind. There are several approaches to solving such equations [7]. In the concrete case of microwave temperature imaging, eq. 7 takes a slightly different form though. This is due to the fact that the brightness temperatures b(f) can practically only be measured at a very restricted number of discrete frequencies, noted f i , i=1, 2, . . . , M, with M around 5 to 7 for the 1D case. This small number of measurements, referred to as sparse data, stands in contrast to the resolution of the temperature profile T(x) one aims to obtain. Therefore, it is practically impossible to discretize T(x) to T(x j ), j=1, 2, . . . , P, where the x j  are the pixel coordinates, in order to find the matrix expression of eq. 1 as b=W·t, and using standard quadrature theory for the solution of eq. 1, such as presented in [8]. This is because in the simplest 1D case, the number P would have already to lie at around 256 or 512 (number of pixels along the 1D profile) in order to reach a practically valuable resolution for the temperature profile T(x). In other words, the matrix formulation would result in a system of only 5 to 7 linear equations with 256 or 512 unknowns T(x j ). Obviously this is a rather hard, if not impossible task to be done reliably, and in the case of a 2D or even 3D reconstruction, this problem would be become even harder. 
     Alternatively, and in contrast to the classical approach to solving Fredholm equations of the first kind, we reformulate eq. 7 in a semi-discrete fashion as follows: 
     
       
         
           
             
               
                 
                   
                     
                       b 
                       
                         i 
                          
                         
                             
                         
                       
                     
                     = 
                     
                       
                         γ 
                         i 
                       
                        
                       
                         
                           ∫ 
                           Ω 
                           
                               
                           
                         
                          
                         
                           
                             
                               
                                 W 
                                 i 
                               
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                             · 
                             
                               T 
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                           
                            
                           
                               
                           
                            
                           
                              
                             x 
                           
                         
                       
                     
                   
                   , 
                   
                       
                   
                    
                   
                     
                       with 
                        
                       
                           
                       
                        
                       i 
                     
                     = 
                     1 
                   
                   , 
                   2 
                   , 
                   … 
                    
                   
                       
                   
                   , 
                   M 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   8 
                 
               
             
           
         
       
     
     where i indexes the M discrete measurement frequencies f i  of the system. This formally quite simple change in problem formulation gives the possibility to solve the described image reconstruction problem for microwave radiometry in an attractive way, in particular for cases where the number M of equations in eq. 8 is small in comparison to the resolution one would like to obtain for the temperature profile T(x). In the following this approach is described in detail. 
     Intuitively, in order to solve eq. 8 for the unknown function T(x), one has to attempt to find a representation for the space of possible temperature profiles T(x), which has to be entirely determined by not much more than the number M of frequencies at which measurements are being carried out, while reflecting the large variability we might find in real-world temperature profiles. In addition, we have to propose algorithms which determine the optimal fit, i.e. the most probable function T(x) for the given measurements b i . The disclosed reconstruction procedure exactly presents a general approach to solve such a problem for microwave radiometry. In general, the disclosed approach to solve eq. 8 for T(x) can be structured into the following two steps:
         1. Construct a function sub-space of L 2  (space of square-integratable functions), denoted Ω D , spanned by an over-complete set of basis functions, D={h k (x)}, k=1, 2, . . . , K, which accounts for the variability of possible temperature profiles T(x). The overcompleteness of the basis functions {h k (x)} means that for any function f(x) from Ω D , L&gt;1 sets {t k,i }, i=1, 2, . . . , L of coefficients exist for which       

     
       
         
           
             
               
                 
                   
                     
                       f 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     ≈ 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                        
                       
                           
                       
                        
                       
                         
                           t 
                           
                             k 
                             , 
                             i 
                           
                         
                         · 
                         
                           
                             h 
                             k 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                       
                     
                   
                   , 
                   
                     i 
                     = 
                     1 
                   
                   , 
                   2 
                   , 
                   … 
                    
                   
                       
                   
                   , 
                   
                     L 
                     . 
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   9 
                 
               
             
           
         
       
         
         
           
              In addition, the set {h k (x)} has to be constructed in a way that for any temperature profile T(x), which can realistically be expected, one set of parameters {t k } out of the L sets of parameters {t k,i }, i=1, 2, . . . , L exists, which approximates the profile T(x) according to eq. 9 with just N&lt;&lt;K non-zero coefficients: 
           
         
       
    
     
       
         
           
             
               
                 
                   
                     
                       T 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     ≈ 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                        
                       
                           
                       
                        
                       
                         
                           t 
                           k 
                         
                         · 
                         
                           
                             h 
                             k 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   10 
                 
               
             
           
         
       
         
         
           
              with the majority of the coefficients t k  being zero, and just about N of the K coefficients having significant magnitude (N&lt;&lt;K). 
             2. Once such a set D has been formed, for any given set of measurements b i  an algorithm has to determine both, the N functions of D and its non-zero weighting coefficients t k  which approximate the temperature profile T(x) optimally by the small weighted sum given in eq. 10. 
           
         
       
    
     The next paragraphs describe these two general steps in more detail. 
     1. Over-Complete Set of Basis Functions for Microwave Radiometry 
     Theoretically, any set of functions of L 2  would be admissible to span the reconstruction subspace of L 2 , noted Ω D , and therefore constitute the basis set, D. Practically though, any prior information about the temperature profile T(x) that can be incorporated into the chosen function set D should be employed in order to be able to implement an efficient and precise reconstruction algorithm. Furthermore, the basis set should be constructed, so that for any real-world temperature profile T(x) just a small subset of N basis functions from the potentially very large set of K basis functions (small in the sense of N&lt;&lt;K) is able to reconstruct T(x). Mathematically this means that real-world temperature profiles T(x) should have sparse representations in the space Ω D . 
     Once the function set D has been constructed, one has still to implement an algorithm which is able to determine the small set of N basis functions from the initial over-complete dictionary of K functions, and its associated non-zero weights t k  which reconstruct an estimate of the unknown temperature profile T(x): 
     
       
         
           
             
               
                 
                   
                     
                       T 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     ≈ 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                        
                       
                           
                       
                        
                       
                         
                           t 
                           k 
                         
                         · 
                         
                           
                             h 
                             k 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   11 
                 
               
             
           
         
       
     
     with just N&lt;&lt;K non-zero coefficients t k . In other words, this algorithms has to determine the set of coefficients {t k } from the sets of coefficients {t k,i }, i=1, 2, . . . , L of equation 9 for which most coefficients are zero. 
     In mathematical terms, we need an algorithm which determines a sparse representation of the temperature profile T(x) characterized by the set of measurements b i  within the subspace Ω D  of L 2 , spanned by the set of functions D. In the next paragraph such algorithms are shown 
     2. Reconstruction Algorithm 
     The initial reconstruction problem of eq. 7 can be rewritten, using the reconstruction approach described in the previous paragraph. Introducing eq. 11 into eq. 8, we get: 
     
       
         
           
             
               
                 
                   
                     b 
                     i 
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       K 
                     
                      
                     
                         
                     
                      
                     
                       
                         t 
                         k 
                       
                       · 
                       
                         γ 
                         i 
                       
                       · 
                       
                         
                           ∫ 
                           Ω 
                           
                               
                           
                         
                          
                         
                           
                             
                               
                                 W 
                                 i 
                               
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                             · 
                             
                               
                                 h 
                                 k 
                               
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                           
                            
                           
                               
                           
                            
                           
                             
                                
                               x 
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   12 
                 
               
             
           
         
       
     
     This problem formulation can also be written in matrix form as: 
         b=A·t,   Eq. 13 
     with 
     
       
         
           
             
               
                 
                   
                     A 
                     
                       i 
                       , 
                       k 
                     
                   
                   = 
                   
                     
                       γ 
                       i 
                     
                     · 
                     
                       
                         ∫ 
                         Ω 
                         
                             
                         
                       
                        
                       
                         
                           
                             
                               W 
                               i 
                             
                              
                             
                               ( 
                               x 
                               ) 
                             
                           
                           · 
                           
                             
                               h 
                               k 
                             
                              
                             
                               ( 
                               x 
                               ) 
                             
                           
                         
                          
                         
                             
                         
                          
                         
                           
                              
                             x 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   14 
                 
               
             
           
         
       
     
     The size of the matrix A is given by the number of measurements M of the system, and by the number of functions K in the over-complete set D. The number of rows is therefore M, while the number of columns is K. 
     Eqs. 12 and 13 allow reinterpreting the temperature image reconstruction in a way which gives direct access to a large range of algorithms perfectly adapted to the problem of finding the optimal weighting coefficients t k  for the profile reconstruction through eq. 10. In fact, eq. 14 can be seen as a transformation F, which transforms the space Ω D  spanned by the basis functions of the set D onto a subspace of R N  (the vector space of N-dimensional vectors of real numbers), denoted Ω R : 
     
       
         
           
             
               
                 
                   
                     
                       F 
                       : 
                       
                         
                           Ω 
                           D 
                         
                         ⋐ 
                         
                           
                             L 
                             2 
                           
                           → 
                           
                             Ω 
                             R 
                           
                         
                         ⊆ 
                         
                           R 
                           N 
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     
                       
                         
                           
                             h 
                             k 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                         → 
                         
                           a 
                           k 
                         
                       
                       = 
                       
                         
                           ∫ 
                           Ω 
                           
                               
                           
                         
                          
                         
                           
                             
                               W 
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                             · 
                             
                               
                                 h 
                                 k 
                               
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                           
                            
                           
                               
                           
                            
                           
                              
                             x 
                           
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     
                       k 
                       = 
                       1 
                     
                     , 
                     2 
                     , 
                     … 
                      
                     
                         
                     
                     , 
                     K 
                     , 
                     
                       
 
                     
                      
                     with 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       W 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     := 
                     
                       ( 
                       
                         
                           
                             
                               
                                 γ 
                                 1 
                               
                               · 
                               
                                 
                                   W 
                                   1 
                                 
                                  
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 γ 
                                 2 
                               
                               · 
                               
                                 
                                   W 
                                   2 
                                 
                                  
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                             
                           
                         
                         
                           
                             KK 
                           
                         
                         
                           
                             
                               
                                 γ 
                                 M 
                               
                               · 
                               
                                 
                                   W 
                                   M 
                                 
                                  
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                             
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   15 
                 
               
             
           
         
       
     
     The resulting subspace Ω R  of R N  can be seen as a vector space spanned by the vectors {a k } which are just the rows in the matrix A of eq. 13. As a last step, we normalize the vectors a k  and therefore the rows of A and reformulate the temperature profile reconstruction of eq. 14 accordingly: 
     
       
         
           
             
               
                 
                   
                     
                       b 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             1 
                           
                           K 
                         
                          
                         
                             
                         
                          
                         
                           
                             t 
                             k 
                             ′ 
                           
                           · 
                           
                             a 
                             k 
                             ′ 
                           
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     with 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         t 
                         k 
                         ′ 
                       
                       = 
                       
                         
                           t 
                           k 
                         
                         · 
                         
                            
                           
                             a 
                             k 
                           
                            
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     and 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       a 
                       k 
                       ′ 
                     
                     = 
                     
                       
                         
                           a 
                           k 
                         
                         
                            
                           
                             a 
                             k 
                           
                            
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   16 
                 
               
             
           
         
       
     
     The symbol ∥a k ∥ stands for the classical vector norm of R N  defined by 
     
       
         
           
             
                 
             
              
             
               
                  
                 
                   a 
                   k 
                 
                  
               
               = 
               
                 
                   
                     
                       ∑ 
                       i 
                     
                      
                     
                         
                     
                      
                     
                       a 
                       i 
                       2 
                     
                   
                 
                 . 
               
             
           
         
       
     
     Eq. 16 is the reformulation of the initial temperature profile reconstruction, characterized by eq. 7, as standard sparse data decomposition using over-complete dictionaries. Solving eq. 17 for a sparse representation aims to determine the few non-zero expansion coefficients t′ k , which result in a good approximation of b. Matching pursuit [3], orthogonal matching pursuit [4], basis pursuit [5], high-resolution matching pursuit [6], etc. are known algorithms for this task. Once this decomposition of b has been determined, the one-to-one correspondence of eq. 15 between the over-complete set of base vectors {a′ k } of Ω R  and the over-complete set of basis functions {h k (x)} of Ω D  enables the temperature profile reconstruction through eq. 10: 
     
       
         
           
             
               
                 
                   
                     T 
                      
                     
                       ( 
                       x 
                       ) 
                     
                   
                   ≈ 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       K 
                     
                      
                     
                         
                     
                      
                     
                       
                         
                           t 
                           k 
                           ′ 
                         
                         
                            
                           
                             a 
                             k 
                           
                            
                         
                       
                       · 
                       
                         
                           
                             h 
                             k 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   17 
                 
               
             
           
         
       
     
     In  FIG. 4 , a functional block-diagram is shown, which describes the different tasks that have to be implemented for the disclosed invention. In  FIG. 5 , some algorithms that can be used in this disclosed invention are shown in relationship to its surrounding functional blocks. The mentioned algorithms are not the only possible algorithms which would solve eq. 16. Common to all of them is that they consist of two parts. First, they have the data part, which considers the actual measurements to determine the optimal temperature profile T(x). And second, there is a sparseness constraint which ensures that just a small number N (small in the sense of N&lt;&lt;K) of the coefficients t′ k  in eq. 17 has significant amplitude (sparse representation) and which, in the case of microwave temperature image reconstruction, plays the role of a regularization term that turns the generally ill-posed inversion problem of microwave temperature image reconstruction into a well-posed inversion problem. The main difference of these “sparseness regularization terms” to the smoothness constraint of Tikhonov regularization (eq. 5), is that they do not over-smooth the final reconstruction result T(x). Let&#39;s shortly describe some of these well-known algorithms. 
     In the notational context of the previous paragraphs, general sparse decomposition using over-complete dictionaries tries to solve the following equation: 
       min|t′∥ 0 , subject to  b=A′·t′,   Eq. 18 
     where ∥.∥ 0  denotes the l 0 -norm. The minimization of the l 0 -norm results in searching for a sparse representation, as the l 0 -norm counts the non-zero expansion coefficients t′ k  in the expansion and tries to minimize their number N under the constrained that the data term is valid. Unfortunately, this problem can only be solved in a combinatorial way, and is therefore computationally very expensive. Therefore one approach consists of replacing the l 0 -norm with a l 1 -norm: 
       min∥t′∥ 1 , subject to  b=A′·t′.   Eq. 19 
     where ∥.∥ 1  denotes the l 1 -norm. The fact of minimizing the l 1 -norm of t′ results in searching for a sparse representation, i.e. just few expansion coefficients t′ k  have significant magnitude. On the contrary to eq. 18, this problem can be solved with linear programming [5]. Interior point methods or simplex methods [13, 14, 15] can be used to solve eq. 19. Matlab&#39;s Optimization Toolbox [12, 15] has standard implementations of these algorithms which are perfectly suited to solve eq. 19. 
     A similar approach as basis pursuit is called basis pursuit denoising. It solves the following equation: 
       min|b−A′·t′∥ 2 +λt′λ 1 ,  Eq. 20 
     where λ is a weighting coefficient, weighting the importance of the regularization term based on the l 1 -norm of the decomposition coefficients t′ k , and ∥.∥ 2  denotes the l 2 -norm. The l 1 -norm component of the minimization problem results in searching for a sparse representation of the initial reconstruction problem. This formulation of sparse data decomposition using over-complete dictionaries can use quadratic programming [16, 17]. Again, this is computationally much more efficient than trying to solve eq. 18, and as for eq. 19, there is a standard implementation of quadratic optimization in Matlab&#39;s Optimization Toolbox [12]. 
     As mentioned, a large number of further known algorithms exist, which can be applied to the presented invention, i.e. to sparse microwave temperature image reconstruction using over-complete dictionaries. A complete overview of these algorithms would go beyond the scope of this text though. Common to all these algorithms is there localization within a concrete system design which would implement the disclosed invention, as was outlined in  FIGS. 4 and 5 . 
     It has to be noted that Thikonov regularization as presented in [2] is not part of the presented class of algorithms, as it does not search for a sparse solution of the inversion problem. This is shown in  FIG. 2 , where almost all coefficients t′ k  have significiant magnitude. In addition, the employed set of Chebychev polynomials does not form a large and over-complete set of basis functions. 
     Detailed Description of Further Embodiments 
     As mentioned, the described reconstruction approach is applicable to 1D, 2D, or 3D reconstruction. The example implementations presented herein are in the context of the reconstruction of 1D temperature profiles for medical applications. Several medical applications could potentially take important profit from an easy-to-use and relatively cheap temperature measuring system. For example, thermal tumor ablation is a medical procedure for cancer treatment whose success is directly related to the exact control of the intra-body temperature distribution. During a thermal tumor ablation treatment, thin needles are introduced percutaneously until their tips reach inside the tumor volume. At the tips, radio-frequency heats the tumor tissue up to a temperature guaranteeing cell death. In this procedure, two characteristics are fundamental for the success of the treatment: On the one hand, all the cancerous tissue has to be killed, and on the other hand, as few healthy tissue as possible should be ablated. As in such a treatment, the killing of the body cells is directly related to reaching or not reaching a particular temperature threshold, the specific system implementation presented herein can significantly improve the security and reliability of the medical treatment plan through real-time monitoring of the internal temperature distribution. 
     Tumor detection is mentioned as a second field of medical application for the presented system implementation. It employs the fact that tumor cells manifest themselves through a small local temperature increase of 1-2° C. Therefore, the temperature monitoring capability of the presented system allows the medical doctor to use temperature as an additional manifestation of cancerous cells during cancer screening. 
     Any concrete implementation of a microwave temperature imaging device has to contain a certain set of functional building blocks. These building blocks were sketched schematically in  FIG. 3 . In the following, two concrete sample implementations of the general functional system are presented. Their main differences lie in the signal processing unit on which the images are reconstructed, and on the concrete reconstruction algorithm being employed for the reconstruction. Concretely, the first sample implementation uses an off-the-shelf personal computer (PC) as the signal processing unit for the image reconstruction, and the employed algorithm is a Matlab implementation of the basis pursuit algorithm as presented in [5]. The second implementation uses an embedded hardware device and the FOCUSS algorithm is used to find the basis pursuit de-noising solution to the reconstruction problem as presented in [18]. 
       FIGS. 6 and 7  illustrate a system comprising a multi-frequency microwave spiral antenna(s) or wave guide antenna(s)  1 , an analogue connection or connections  6 , a Dicke-null-balancing radiometer  7 , a PC data bus, such as USB, PCI, FireWire, or others  8  and a PC. The system according to  FIG. 6  comprises an off-the-shelf-PC  5 ,  9  with Matlab drivers for the analogue-to-digital data converter (Iotech DaqBoard/2000 Series)  15  and with the reconstruction algorithm implemented in Matlab. The system according to  FIG. 7  comprises a PC  5 ,  9 ,  15  with Iotech DaqBoard/2000 Series, corresponding Matlab drivers and reconstruction algorithms implemented in Matlab. In  FIGS. 6 and 7 , the building blocks of the PC-based sample implementation are schematically shown, while the embedded sample implementation is drawn in  FIGS. 11 and 12 . Before presenting the implementation specific aspects of the devices, the common parts of both implementations are presented. 
     Both systems comprise one single (1D temperature profile) multi-frequency spiral microwave antenna  1  as disclosed by Jacobsen and Stauffer in [9]. The employed antenna operates at 7 frequencies in the range of 0.5 to 3.75 GHz. The signal sensed by the antenna passes through an analogue connection  6  to a Dicke null-balancing radiometer  7 , as described by Jacobsen and al. in [10]. The Dicke radiometer detects sequentially the signal contributions at the different frequencies of the multi-frequency antenna of the system. The resulting sequential analogue signals from the radiometer are directly related to the brightness temperatures at the corresponding frequencies, and these brightness temperatures are directly related to the real intra-body temperature distribution through eq. 8. In order to reconstruct a one dimensional profile of intra-body temperatures from the brightness temperatures, the reconstruction procedure disclosed herein is getting applied to the brightness temperatures, i.e. to the outputs from the analogue-to-digital converter  3 ,  10 ,  15 . 
     As mentioned in the general description, any specific reconstruction implementation needs to specify two aspects: first, the over-complete set of basis functions {h k (x)} has to be specified, and second, the optimization algorithm has to be chosen. It is important to note that in the present context of 1D temperature profile reconstruction, the spatial variable x parameterizes just a 1D space Ω, i.e. is represented by one-component vectors. The size of the imaging space Ω is given by the maximal sensing depth d max  of the microwave antenna, such as e.g. 10 cm for the cases presented herein. 
     In real-world medical applications for microwave temperature monitoring, such as cancer detection or thermal tumor ablation monitoring, bell-shaped temperature distributions can be expected. Therefore, we propose to use 1D B-splines at different scales and positions to build up the over-complete set D of K basis functions. 
     Herein, B-splines of order N d =3 (cubic B-splines) and at N s =30 consecutive scales are used to build the over-complete set D of base functions. The overcompleteness of D arises from the different scales of the B-splines. At the lowest scale, 4 B-splines cover the whole reconstruction space Ω, while at each higher scale, one spline is added. Therefore, at each scale s, N p(s) =s+4 B-splines at different positions are added to the set D, and the whole over-complete set of basis functions can be parameterized as: 
         D={h   k ( x )}={β s,p ( s ) N     d   ( x )},  Eq. 21 
     N d =3, 
     s=0, 1, . . . , 29,
 
p(s)=1, 2, . . . , s+4.
 
with
 
     
       
         
           
             
               
                 
                   
                     
                       β 
                       
                         s 
                         , 
                         
                           p 
                            
                           
                             ( 
                             s 
                             ) 
                           
                         
                       
                       3 
                     
                     ( 
                     
                         
                     
                      
                     x 
                     ) 
                   
                   = 
                   
                       
                   
                    
                   
                     { 
                     
                         
                     
                      
                     
                       
                         
                           
                             
                               
                                 2 
                                 3 
                               
                               - 
                               
                                 
                                    
                                   x 
                                    
                                 
                                 2 
                               
                               + 
                               
                                 
                                   
                                      
                                     x 
                                      
                                   
                                   3 
                                 
                                 2 
                               
                             
                             , 
                           
                         
                         
                           
                             
                               
                                 if 
                                  
                                 
                                     
                                 
                                  
                                 0 
                               
                               ≤ 
                               
                                  
                                 
                                   
                                     
                                       
                                         
                                           ( 
                                           
                                             x 
                                             - 
                                             
                                               
                                                 d 
                                                 max 
                                               
                                               · 
                                               
                                                 
                                                   ( 
                                                   
                                                     
                                                       p 
                                                        
                                                       
                                                         ( 
                                                         s 
                                                         ) 
                                                       
                                                     
                                                     - 
                                                     2 
                                                   
                                                   ) 
                                                 
                                                 
                                                   
                                                     N 
                                                     
                                                       p 
                                                        
                                                       
                                                         ( 
                                                         s 
                                                         ) 
                                                       
                                                     
                                                   
                                                   - 
                                                   3 
                                                 
                                               
                                             
                                           
                                           ) 
                                         
                                         · 
                                       
                                     
                                   
                                   
                                     
                                       
                                         
                                           
                                             N 
                                             
                                               p 
                                                
                                               
                                                 ( 
                                                 s 
                                                 ) 
                                               
                                             
                                           
                                           - 
                                           3 
                                         
                                         
                                           d 
                                           max 
                                         
                                       
                                     
                                   
                                 
                                  
                               
                               ≤ 
                               1 
                             
                             , 
                           
                         
                       
                       
                         
                           
                             
                               
                                 
                                   ( 
                                   
                                     2 
                                     - 
                                     
                                        
                                       x 
                                        
                                     
                                   
                                   ) 
                                 
                                 3 
                               
                               6 
                             
                             , 
                           
                         
                         
                           
                             
                               
                                 if 
                                  
                                 
                                     
                                 
                                  
                                 1 
                               
                               ≤ 
                               
                                  
                                 
                                   
                                     
                                       
                                         
                                           ( 
                                           
                                             x 
                                             - 
                                             
                                               
                                                 d 
                                                 max 
                                               
                                               · 
                                               
                                                 
                                                   ( 
                                                   
                                                     
                                                       p 
                                                        
                                                       
                                                         ( 
                                                         s 
                                                         ) 
                                                       
                                                     
                                                     - 
                                                     2 
                                                   
                                                   ) 
                                                 
                                                 
                                                   
                                                     N 
                                                     
                                                       p 
                                                        
                                                       
                                                         ( 
                                                         s 
                                                         ) 
                                                       
                                                     
                                                   
                                                   - 
                                                   3 
                                                 
                                               
                                             
                                           
                                           ) 
                                         
                                         · 
                                       
                                     
                                   
                                   
                                     
                                       
                                         
                                           
                                             N 
                                             
                                               p 
                                                
                                               
                                                 ( 
                                                 s 
                                                 ) 
                                               
                                             
                                           
                                           - 
                                           3 
                                         
                                         
                                           d 
                                           max 
                                         
                                       
                                     
                                   
                                 
                                  
                               
                               ≤ 
                               2 
                             
                             , 
                           
                         
                       
                       
                         
                           
                             0 
                             , 
                           
                         
                         
                           
                             otherwise 
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   22 
                 
               
             
           
         
       
     
     In  FIG. 8 , the cubic B-spline basis functions for scale 0 and 9 are drawn in dependence on the profile depth in cm. The cardinality of D, |D|=K, is given by: 
     
       
         
           
             
               
                 
                   K 
                   = 
                   
                     
                        
                       D 
                        
                     
                     = 
                     
                       
                         
                           ∑ 
                           
                             s 
                             = 
                             0 
                           
                           
                             
                               N 
                               s 
                             
                             - 
                             1 
                           
                         
                          
                         
                           N 
                           
                             p 
                              
                             
                               ( 
                               s 
                               ) 
                             
                           
                         
                       
                       = 
                       
                         
                           
                             ∑ 
                             
                               s 
                               = 
                               0 
                             
                             
                               
                                 N 
                                 s 
                               
                               - 
                               1 
                             
                           
                            
                           s 
                         
                         + 
                         4 
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   23 
                 
               
             
           
         
       
     
     Therefore in the present case, we have |D|K=555. In accordance with the disclosed reconstruction procedure, we have to transform the base functions of D by the transformation given in eq. 15. In order to be able to do this, we need to specify the explicit weighting functions characterizing the used antenna. In the case of 1D reconstruction, an analytic model of the weighting function is given in [2]: 
     
       
         
           
             
               
                 
                   
                     
                       
                         W 
                         i 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         1 
                         
                           d 
                           i 
                         
                       
                       · 
                       
                          
                         
                           
                             - 
                             x 
                           
                           / 
                           
                             d 
                             i 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   24 
                 
               
             
           
         
       
     
     where the parameters d i  are called the antenna sensing depths at frequencies f i . For the employed spiral antenna and the chosen frequencies, the sensing depths lie at {0.9 cm, 1.2 cm, 1.5 cm, 1.8 cm, 2.1 cm, 2.4 cm, and 2.7 cm}. For this weighting function and for the chosen B-spline basis functions, an analytical solution for the transformation integral of eq. 15 exists, which can be determined easily using the following recursive formula: 
     
       
         
           
             
               
                 
                   
                     ∫ 
                     
                       
                         
                           x 
                           n 
                         
                         · 
                         
                            
                           ax 
                         
                       
                        
                       
                          
                         x 
                       
                     
                   
                   = 
                   
                     
                       
                         1 
                         a 
                       
                       · 
                       
                         x 
                         n 
                       
                       · 
                       
                          
                         ax 
                       
                     
                     - 
                     
                       
                         n 
                         a 
                       
                        
                       
                         ∫ 
                         
                           
                             
                               x 
                               
                                 n 
                                 - 
                                 1 
                               
                             
                             · 
                             
                                
                               ax 
                             
                           
                            
                           
                             
                                
                               x 
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   25 
                 
               
             
           
         
       
     
     This equation can be verified easily using partial integration, or it can be found in mathematical reference books such as in [23]. Therefore the vectors {a k ′} can be calculated analytically. If the transformation integral of eq. 15 cannot be solved analytically because the weighting functions W i (x) have a more complex form than given in eq. 24, or because electromagnetic simulations have been used to obtain a discrete model of these weighting functions, numerical integration has to be employed. 
     The reconstruction process can formally be summarized for the sample implementations presented in this section as follows: 
         b=A′·t′,   Eq. 26 
     where b is the vector whose components are the measured brightness temperatures at the 7 frequencies. The matrix A′ is constituted of K=555 normalized vector columns with 7 components each, noted a k ′, and t′ are the K=555 expansion coefficients one likes to recover. Therefore, eq. 14 becomes explicitely: 
     
       
         
           
             
               
                 
                   
                     
                       
                         A 
                         
                           i 
                           , 
                           k 
                         
                         ′ 
                       
                       = 
                       
                         
                           A 
                           
                             i 
                             , 
                             k 
                           
                         
                         
                            
                           
                             a 
                             k 
                           
                            
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     with 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         A 
                         
                           i 
                           , 
                           k 
                         
                       
                       = 
                       
                         
                           ∫ 
                           Ω 
                         
                          
                         
                           
                             
                               
                                 γ 
                                 i 
                               
                               
                                 d 
                                 i 
                               
                             
                             · 
                             
                               exp 
                                
                               
                                 ( 
                                 
                                   x 
                                   
                                     d 
                                     i 
                                   
                                 
                                 ) 
                               
                             
                             · 
                             
                               
                                 β 
                                 k 
                                 3 
                               
                                
                               
                                 ( 
                                 x 
                                 ) 
                               
                             
                           
                            
                           
                              
                             x 
                           
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     
                       i 
                       = 
                       0 
                     
                     , 
                     1 
                     , 
                     2 
                     , 
                     … 
                      
                     
                         
                     
                     , 
                     
                       
                         6 
                          
                         
                             
                         
                          
                         and 
                          
                         
                             
                         
                          
                         k 
                       
                       = 
                       0 
                     
                     , 
                     1 
                     , 
                     2 
                     , 
                     … 
                      
                     
                         
                     
                     , 
                     554. 
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   27 
                 
               
             
           
         
       
     
     Two concrete algorithms to solve eq. 26 for a sparse solution (just N&lt;&lt;K components of t′ are non-zero) of the unknowns t′ are presented in the context of two concrete sample implementations in the following two sections. Once the parameters t′ have been determined, the temperature profile T(x) is calculated as 
     
       
         
           
             
               
                 
                   
                     T 
                      
                     
                       ( 
                       x 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         0 
                       
                       554 
                     
                      
                     
                       
                         
                           t 
                           k 
                           ′ 
                         
                         
                            
                           
                             a 
                             k 
                           
                            
                         
                       
                       · 
                       
                         
                           
                             β 
                             k 
                             3 
                           
                            
                           
                             ( 
                             x 
                             ) 
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   Eq 
                   . 
                   
                       
                   
                    
                   28 
                 
               
             
           
         
       
     
     The resulting reconstructed temperature profile T(x) is finally visualized on a computer screen  5  for visual inspection by the medical doctor or industrial investigator. 
     In the next two sections, the implementation specific aspects of the presented sample implementations will be presented. 
     1. PC-based Sample Implementation 
     The sample implementation outlined in  FIGS. 6 and 7  is based on an off-the-shelf PC  9 . It also comprises a PCI analogue-to-digital converter  15 , such as Iotech&#39;s DaqBoard\2000 [11] with the provided Matlab [12] drivers, plugged in a PCI-slot  8  of the core PC. This converter converts the analogue signals from the Dicke null-balancing radiometer  7 , which represents the voltage encoded brightness temperatures, b, into their digital representations which are directly available inside the Matlab implementation of the reconstruction algorithm. This algorithm implements the general building blocks outlined in  FIG. 5 , by calculating the matrix A′ through eq. 27, and determining the expansion coefficients t′ by basis pursuit [5]. Basis pursuit solves in the present context the following equation through the interior point algorithm already implemented in Matlab [14, 15]: 
       min∥t′∥ 1 , subject to  b=A′·t′.   Eq. 29 
     The resulting solution vector t′ represents the expansion coefficients to be used for the temperature profile estimation by eq. 28. 
     In  FIG. 9 , three sample results are shown for the presented sample implementation, which used 100 iterations, and a convergence precision of 10 −3  for the interior point algorithm of the Matlab Optimization Toolbox [12]. The reconstruction results are represented in broken lines and the references in full lines. We see the clear advantage of the disclosed reconstruction approach compared to current state-of-the-art microwave temperature profile reconstruction shown in  FIG. 2 . 
     2. Embedded Sample Implementation 
       FIG. 10  shows an embedded system for microwave image reconstruction which is connected to a multifrequency spiral antenna  1  via an analogue connection  6  and comprises a dicke null-balancing radiometer  7 , an analogue-to-digital converter from Orsys, such as ORS-116, and the embedded development board C6713 Compact from Orsys connected via a firewire connection  12  to a core PC with Unibrain firewire drivers, wherein via said connection  12  images are sent to the core PC and configuration to the embedded system. 
     The system according to  FIG. 11  comprises a PC with firewire drivers and C++ API from Unibrain  5 ,  13 , a microwave spiral antenna  1  and an embedded system  7 ,  10 ,  11  for microwave image reconstruction including a Dicke null-balancing radiometer, an analogue-to digital converter from Orsys and a C6713 Compact development board from Orsys. The embedded system  7 ,  10 ,  11  is connected to the PC  5 ,  13  via a firewire connection  12  and via an analogue connection  6  to the antenna  1 . 
     For the embedded sample implementation as outlined in  FIGS. 10 and 11 , the reconstruction algorithm is implemented on an embedded system  11  from Orsys [21]. More precisely, the radiometer  7  outputs, representing the voltage-encoded brightness temperatures sensed by the multi-frequency spiral antenna, are digitized by the ORS-116 data converter  10  from Orsys [21]. The resulting vector, b, with the brightness temperature components is being employed to reconstruct the temperature profiles on the C6713 Compact  11  from Orsys [21]. The reconstruction result T(x) is transmitted through the FireWire port  12  of the C6713 Compact device to the host computer  13 . On the host PC  13 , the FireWire drivers by Unibrain [22] are being employed to receive the temperature profiles from the C6713 Compact device  11 . Finally, these temperature profiles are getting visualized on the computer screen  5  using the software packages vtk from Kitware [19] and Qt from TroUtech [20]. 
     The reconstruction algorithm implemented on the embedded C6713 Compact device  11  from Orsys [21], calculates the matrix A′ through the mathematical expression given in eq. 27. The FOCUSS algorithm detailed in [18] is used to determine the expansion coefficients t′ for the reconstruction of the temperature profile T(x) through eq. 28. This algorithm solves the following expression iteratively as described in  FIG. 12 : 
       min∥b−A′·t′∥ 2 +λ∥t′∥ 1 .  Eq. 30 
     Eq. 30 represents the so-called “basis pursuit de-noising” approach to data decomposition using over-complete dictionaries. In  FIG. 13 , the reconstruction results for this sample implementation are shown for 500 iterations, λ=0.000001, and where the initial estimate of t′ 0  (before any iteration) is a null vector. The reconstruction results are represented in broken lines and the references in full lines. The resulting coefficients t′ determine the temperature profile T(x) by the weighted sum of cubic B-splines as given in eq. 28. 
     REFERENCES 
     
         
         [1] B. Stec, A. Dobrowolski, and W. Susek: “Multifrequency microwave thermography for biomedical applications” IEEE Transactions on Biomedical Engineering, Vol 51(3), March 2004 
         [2] S. Jacobsen and P. R. Stauffer: “Nonparametric 1-D temperature restoration in lossy media using Tikhonov regularization on sparse radiometry data”, IEEE Transactions on Biomedical Engineering, Vol. 50(2), February 2003 
         [3] S. Mallat and Z. Zhang: “Matching pursuit with time-frequency dictionaries”, IEEE Transactions on Signal Processing, Vol. 41, December 1993 
         [4] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad: “Orthogonal matching pursuits: recursive function approximation with applications to wavelet decomposition”, Proceedings of the 27 th  Asilomar Conference in Signals, Systems, and Computers, 1993 
         [5] S. S. Chen, D. L. Donoho, and M. A. Saunders: “Atomic decomposition by basis pursuit”, SIAM J. Sci. Comput., Vol 20(1), 1998 
         [6] S. Jaggi, W. C. Karl, St. Mallat, and A. S. Willsky: “High resolution pursuit for feature extraction”, Technical Report, MIT, November 1996 
         [7] L. M. Delves, and J. L. Mohamed: “Computational methods for integral equations”, Cambridge University Press, 1985 
         [8] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery: “Numerical Recipes in C”, Cambridge University Press, 1992 (Second Edition) 
         [9] S. Jacobsen and P. R. Stauffer “Multifrequency Radiometric Determination of Temperature Profiles in a Lossy Homogeneous Phantom Using a Dual-Mode Antenna with Integral Water Bolus”, IEEE Transactions on Microwave Theory and Techniques, vol. 50(7), pp. 1737-1746, July 2002. 
         [10] S. Jacobsen, P. R. Stauffer, and D. Neuman: “Dual-mode Antenna Design for Microwave Heating and Noninvasive Thermometry of Superficial Tissue Disease”, IEEE Transactions on Biomedical Engineering, vol. 47, pp. 1500-1509, November 2000. 
         [11] http:/www.iotechcom 
         [12] http://www.matlab.com 
         [13] G. B. Dantzig, A. Orden, and P. Wolfe: “Generalized Simplex Method for Minimizing a Linear from Under Linear Inequality Constraints”, Pacific Journal Math, vol. 5, pp. 183-195. 
         [14] S. Mehrotra: “On the implementation of a Primal-Dual Interior Point Method”, SIAM Journal on Optimization, vol. 2, pp. 575-601, 1992. 
         [15] Y. Zhang: “Solving Large-Scale Linear Programs by Interior-Point Methods under the MATLAB Environment”, Technical Report TR96-01, Department of Mathematics and Statistics, University of Maryland, Baltimore County, Baltimore, Md., July 1995 
         [16] T. F. Coleman, and Y. Li: “A Reflective Newton Method for Minimizing a Quadratic Function Subject to Bounds on some of the Variables”, SIAM Journal on Optimization, vol. 6(4), pp. 1040-1058, 1996. 
         [17] P. E. Gill, W. Murray, and M. H. Wright: “Practical Optimization” Academic Press, London, UK, 1981. 
         [18] B. D. Rao, K Engan, S. F. Cotter, J. Palmer, and K Kreutz-Delgado: “Subset Selection in Noise Based on Diversity Measure Minimization”, IEEE Transactions on Signal Processing, vol 51(3), pp. 760-770, March 2003 
         [19] http://www.vtk.org 
         [20] http://www.trolltech.com 
         [21] http://www.orsys.de 
         [22] http://www.unibrain.com 
         [23] I. N. Bronstrein, K. A. Semendjajew, G. Musiol, and H. Mühlig “Taschenbuch der Mathematik”
       Verlag Harri Deutsch, 2 nd  edition, 1995